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December,  1966 

Chairman:   Dr.  Andrew  P.  Sage 

Major  Department:   Nuclear  Engineering  Sciences 

The  main  objective  of  this  study  is  to  investigate 
optimum  control  of  the  dynamics  of  nuclear  systems  and 
associated  computational  problems.   Two  systems  are  con- 
sidered:  1)  Nuclear  Reactor  System,  2)  Nuclear  Rocket 
Engine  System.   For  the  first  system,  three  problems  are 
analyzed:   a)  the  neutron  density  or  reactor  power  is  trans= 
ferred  from  one  steady  state  condition  to  another  with  a 
quadratic  constraint  on  control  rod  movement.   By  varying 
the  operating  time  and  placing  an  inequality  constraint  on 
reactivity  a  form  cf  minimum  time  control  is  obtained,  b) 
integral  of  quadratic  form  of  error  in  desired  neutron 
density  is  minimized  with  quadratic  integral  constraint 
on  control  rod  movement,  c)  a  fixed  configuration  optimum 
closed-loop  control  system  is  realized  by  letting  control 
rod  nrcvement  be  a  proportional  plus  integral  function 
of  neutron  density  error. 
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The  two-point  boundary  value  problems  (TPBVP) ,  re- 
sulting by  the  application  of  optimization  theory  to  point 
reactor  models  are  obtained.   A  modified  gradient  technique 
is  used  to  exploit  the  hybrid  capabilities  of  modern  analog 
computers.   The  results  are  compared  with  the  solutions 
of  this  problem  by  the  numerical  technique  of  quasilinear= 
ization.   In  all  of  these,  the  convergence  is  obtained  in 
no  more  than  four  iterations. 

In  the  examples  considered  here,  the  cost  function 
is  higher  for  the  closed  loop  case  as  opposed  to  the  open 
loop  case.   This  is  in  general  agreement  with  physical 
insight  and  previous  results,  since  an  assumed  form 
of  closed  loop  controller  is  used.   The  coefficients  A 
and  B  in  the  control  law  are  highly  dependent  upon  the  final 
time  tf,  the  system  parameters,  and  the  initial  and  desired 
final  state.   Variation  in  system  parameters  may,  however, 
render  the  closed  loop  control  better  than  the  open  loop 
control  which  does  not  change  with  changing  system  parameter; 

For  the  case  of  the  nuclear  rocket  engine  system, 
only  a  quadratic  constraint  on  total  reactivity  is  con= 
sidered.   The  nuclear  rocket  engine  is  similar  to  the  solid- 
core  25,000  megawatt,  million  pound  thrust  model  developed 
by  Smith  and  Stenning.   Here,  the  results  obtained  for 
various  start-up  times  and  various  desired  power  levels  in 
terms  of  control  needed  and  required  cost,  appear  very 
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reasonable.   Also  it  was  observed  that  if  the  pressure  and 
the  temperature  tirre  constants  were  increased  frotr  their 
normal  values ,  the  temperature  response  became  somewhat 
sluggish.   Possible  improvements  to  this  situation  are 
suggested. 
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CHAPTER  I 

INTRODUCTION 

Just  as  all  other  fields  of  scientific  thought  have 
undergone  a  major  change  in  recent  years,  so  also  has  the 
field  of  rocketery.   From  the  Chinese  rockets  of  1232 
A.D.  to  those  of  the  Germans  during  the  second  world  war, 
man  has  made  a  "case"  for  nuclear  rocket  propulsion  as  a 
vital  factor  in  the  success  of  future  space  travels   This 
has  been  pointed  out  in  the  recent  literature  (1,2,3,4) 
and  accepted  in  the  field  today.-   Its  acceptance  is  fur= 
ther  strengthened  by  the  successful  testing  of  the  KIWI 
and  NERVA  reactors. 

This  sudden  interest  in  space  travel  has  led  to 
thoughts  about  improved  system  performances,  as  a  result 
of  which  a  very  instructive  book  on  nuclear  rocket  system 
optimization  by  Bussard  (5)  has  appeared.   An  outgrowth 
of  this  interest  is  the  optimal  control  problem  wherein 
the  control  is  such  as  to  give  the  "best  performance". 
From  the  nuclear  literature,  however,  it  is  clear  that 
the  subject  of  optimal  control  of  nuclear  systems  has 
received  little  interests   Recently  a  num.ber  of  papers 
have  appeared  in  this  area,  ranging  from  optimal  reactor 
shutdown  programs  for  control  of  xenon  poisoning  (6,7) , 


and  synthesis  of  optimal  nuclear  reactor  systems  (8,9)  to 
optimal  control  of  nuclear  reactor  systems  (10. 11 ) . 
More  recently  Mohler  and  Weaver  have  reported  some  work 
in  the  application  of  modern  optimal  control  theory  to 
the  dynamics  of  nuclear  rocket  systems  (ll.,12)  but  much 
still  remains  to  be  done. 

The  mathematical  theory  used  for  the  study  of 
optimization  problems  is  the  calculus  of  variations.   It 
is  of  interest  to  note  that  Goddard  (l_3)  recognized  this 
technique  as  an  important  tool  in  the  performance  analysis 
of  rockets  in  his  early  paper  published  40  years  ago. 
There  are  many  possible  variational  methods  fcr  minimizing 
or  maximizing  a  functional  over  a  function  space.   Among 
the  commonly  used  methods  in  control  system  design  are: 

(1)  calculus  of  variations 

(2)  maximum  principle 

(3)  dynamic  programming 

In  all  cases  the  object  is  to  find  the  optimum  con- 
trol law  such  that  the  given  functional  of  the  performance 
indices  is  minimized  or  maximized.   These  methods  are  re= 
lated  (5)  in  that  the  rr.aximum  principle  can  be  derived  from 
the  calculus  of  variations   as  well  as  from  the  method 
of  dynamic  programming^    Pontryagin's  maximum 


principle  is  a  very  elegant  method,  which  is  applicable 
to  a  broader  class  of  problems  than  the  calculus  of 
variations.   An  outcome  of  this  principle  is  the  so  called 
two»point  boundary  value  problem,  abbreviated  as  "TPBVP',' 
with  split  boundary  conditions. 

This  method  will  be  applied  in  the  study  of  dynamics 
of  nuclear  systems.   Since  the  dynamic  response  of  nuclear 
systems  is  represented  by  a  set  of  nonlinear,  time-varying 
differential  equations,  conventional  analysis  is  usually 
too  difficult  to  give  a  closed  form  solution.   Develop- 
ment of  efficient  computer  techniques  is  therefore  necessary. 

A  numerical  approach  which  is  currently  quite 
popular  is  the  method  of  gradients  (5).   It  proceeds  by 
solving  a  sequence  of  sub-optimal  problems  with  the  prop- 
erty that  each  successive  set  of  solution  functions  yields 
an  improved  value  for  the  functional  being  optimized. 

Another  numerical  technique  of  considerable  interest 
is  the  approximating  process  of  quasilinearization  which 
produces  a  sequence  of  approximations  x^,  n  =  l,2,..c  that 
converge  quadratically  (14)  to  the  solution  of  the  ori= 
ginal  two=point  boundary  value  problem.   This  technique 
was  initially  suggested  by  Bellman  and  Kalaba  (15)  as  a 
modification  of  the  Newton-Raphson-Kantorovich  method. 
This  converts  the  problem  of  solving  a  two-point  boundary 
non-linear  differential  equation  into  a  problem  of  solving 
a  sequence  of  two-point  boundary  linear  differential  equation 


problems.   The  solution  of  the  latter  equation  is  con- 
siderably easier  to  find. 

For  the  study  of  optimization  problems,  a  precise 
formulation  of  the  performance  index  is  necessary.   The 
performance  index  is  assumed  to  be  given  by: 

ftf 


J(Xo,tf,u) 


0  (x,u,t)dt 


where  0  is  a  scalar,  a  function  of  the  state  vector  x,  the 
control  vector  u  and  time. 

Corresponding  to  this  performance  index,  a  control  law 
is  found  which  is  the  best  in  the  sense  of  maximiizing 
the  performance  index. 

The  performance  indices  considered  in  this  study  are: 

■f     ( 


i) 


J  =  1/2 


t: 


dt 


a)  n(tf)  =  n^i 

b)  T(tf)  =  Td 

This  performance  criterion  may  be  interpreted  as  an  optimi- 
zation of  control  rod  movement,  limiting  erratic  flux  changes 
(10).   Clearly  it  tends  to  limit  the  maximum  reactor  period. 

J  =  1/2   (  ^     ^c2  +  a  (n=nd)^    dt 
0     L 

Integral  of  the  error  squared  of  the  reactor  power  and  its 
desired  value  is  minimized  along  with  the  integral  of  the 
square  of  the  reactivity  function,  a  is  the  error  weight- 
ing term 
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iii)  J  =  1/2   j  ^   [  ^c^  +  a(n-nd)^  ]  dt 


where   ^c  =  a  +  ACn-n^) 


a   and  a  are  constants  to  be  determined, 
iv)  J  =  1/2    )  ^  [  ^c^  +  a(n-n.)^J  dt 


where     P^   =  Atn-n^j)  +    b)    (n-n^)  dt 

o 


A  and  B  are  constants  to  be  determined. 

The  last  two  criteria  give  closed  loop  control  and  since 
a  specific  form  is  chosen  for  the  controller,  it  will 
henceforth  be  referred  to  as  the  specific  optimal  control 
problem. 

The  above  four  controls  are  analyzed  and  compared 
with  one  another  for  two  reactor  models.  The  two  models 
analyzed  are: 

(1)  Point  model  of  a  nuclear  reactor  system. 

(2)  Point  model  of  a  nuclear  rocket  engine. 

Chapter  II  contains  the  description  of  the  first  model. 
Chapter  III  concerns  the  second  model  and  follows  along 


the  same  pattern  as  Chapter  II.   Chapter  IV  includes 
a  discussion  of  results  and  suggestions  for  further 
study. 


CHAPTER  II 
OPTIMAL  COKTRCL  CF  A  NUCLEAR  REACTOR  SYSTEM 
A.    Introduction  and  Model  Description 

Need  frequently  arises  in  nuclear  reactors  as  well 
as  in  nuclear  power  plants  for  power  level  changes.   These 
changes  are  generally  accomplished  by  semi-automatic 
means  in  small  research  reactors.   In  large  power  re- 
actors and  in  nuclear  rocket  engines,  however,  automatic 
control  is  desired.   If  this  control  is  optimal  in  some 
specified  sense,  economic  or  other  advantages  could  be 
gained. 

The  reactor  model  considered  here  is  that  of  a  one 
delayed-neutron  group  without  temperature  feedback  and 
without  spatial  and  angular  dependence,  described  by  the 
following  system  of  differential  equations: 

n  =  -^  (P  -  /3  )n  +  /.c  (2.1) 

c  =  -f-  n  -  AC  (2.2) 

where  n  and  c,  the  neutron  density  and  precursor  popula- 
tion are  the  state  variables  and  the  reactivity  ^j      is 


the  control  variable.   The  dot  represents  derivative 
with  respect  to  time. 

A  ,  is  the  average  neutron  lifetime 

/.  ,  is  the  decay  constant,  and 

/3  ,    is  the  delayed  neutron  function. 

The  system  is  at  steady  state  for  t  <  o  and  the 
initial  conditions  are:   n(o)  =  n^  and  c(o)  =  c^. 
'     It  is  desired  to  optimize  system  performance  such  as  to 
minimize  the  following  integrals. 


1)   J  =  1/2  f        Pc^it)      dt 
•^0 

with  n(tf)  =  n^  and   I  P^  I    <  Pmax 


2)   J  =  1/2  [^    [Pc^   +  o{n-n^f] 


dt 


in  fixed  time  tj 


3)   J  =  1/2  J^    \^     Pc2  +  a(n-nd)2  1 


dt 


where   p^  =  ^  +  A(n=n^) 
in  fixed  time  t£ 


4)   J  =  1/2  J^   [  Pc2  +  a  (n-n^)^  ]  dt 


where  p^  -   A{n-nci)  +  B  f       (n-n^)   dt 


in  fixed  time  tf. 


1 .   Staterrent  of  the  K^axitnum  Principle 

Consider  the  control  processes  to  be  described  by 
a  vector  differential  equation: 

X  =  f (x,u,t)   ,   x(c)  =  Xo  (2.3) 

where  f(x,u,t)  is  a  vector  with  coordinates  f^(x,u,t), 

f2(x.u,t) fn(x.U,t). 

The  control  functions  have  to  be  piecewise  continuous 
and  can  be  restricted  by  the  inequalities: 

gj(u^,...u^)  <  0    (j  =  l....m)  (2.4) 


Assume  a  perforrrance  criterion  of  the  type: 

tf 
J=  0[x(tf).u(tf),  ^^]+  /  <P(x,Z,t)      dt       (2.5) 

where  "J"  is  generally  referred  to  as  the  performance 
criterion  or  performance  index  or  cost  function. 

The  Hamiltonian  of  the  system  is  now  defined  by. 

H(x,u,p.t)  =  0(x,u,t)  +  P^(t)  f(x,u.t)        (2.6) 


This  Hamiltonian  becomes  the  Hamiltonian  of  classical 
mechanics  only  if  the  control  u  becomes  optimal  control 
u*  defined  by  equation  (2.9). 
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The  variable  appearing  in  Equation  (2.6)  is  called 
the  auxiliary  vector  or  the  adjoint  vector  or  momentum 
function.   The  system  equations  can  now  be  written  in  the 
following  manner. 


aH(x.u.p.t) 

a  Pi 


(2.7) 


Auxiliary  or  adjoint,  or  costate  functionsp^^  are  defined  by 

(2.8) 


p.(t)  =  -  '^H(x.u,p,t) 
d  ^x 


The  Hamiltonian  is  now  minimized  with  respect  to  the  control 
variable  u,  such  that 


HO(x.p,t)  =  Min  H(x,p,u,t)  (2.9) 

where  u°(x,p,t)  minimizes  H  and  the  symbol^x  refers  to 

the  admissible  control  space.   If  rx  is  infinite  (2.9)  be- 
comes  -^  =  0.   The  boundary  conditions  are  given  by  : 


x(to)  =  Xq 


(2.10) 


and  P(t^)  =   V   a  (x.u.t) 
2.   Specific  Optimal  Control 


where  (  Vx  ^  )^  = 
t  =  tr 


be 


The  desirability  of  closed  loop  control  laws  has  led 
to  the  development  of  "specific  optimal  control." 
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The  specific  optimal  control  (S.O.C.)  problem  is 
defined  in  the  following  manner  (17)  :  given  a  plant  with 
a  state  equation  of  the  form 

X  =  T{x,u,t)  ,  x(to)  =  Xq  (2.11) 

It  is  desired  to  determine  the  unknown  parameters  in  a 
control  law  of  the  form 

u  =  h{y,b)  .    b  =  0  (2.12) 

In  the  above  equation  y  is  a  p-dimensional  vector  which  is 
a  known  function  of  the  state  and  F  is  a  q-dimensional 
constant  vector  to  be  determined  such  that  an  index  of 
performance  of  the  form 

^f   - 
J(u)  =    /•   g(x,u)  dt  (2.13) 

Jo 

is  minimized. 

The  specific  optimal  control  problem  can  be  solved 
by  a  simple  re-formulation  of  the  above  problem,,   The  re- 
formulation is  carried  out  by  substituting  equation  (2.12) 
into  equations  (2.11)  and  (2.13),  leading  to 

x  =  f (x,h(y,b),t) 

tf    _  ■  _  _ 
J(b)  =   r    g(x.h(y.b))  dt 
•'o 
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Since  y  is  a  known  function  of  x  and  b  is  a  constant  vector 


X  =  f (x.b.t) 

^f 
J(b)  =    f       g(x.b)dt  (2.14) 

•'o 


and  b  =  0 

Thus  the  S.O.C.  problem  has  been  embedded  in  the  pro- 
blem of  finding  the  vector  b,  which  minimizes  the  cost 
function  and  differential  constraints  imposed  by  equation 
(2.14).   The  above  is  now  in  the  form  of  the  typical  opti- 
mal control  problem.   Using  Pontryagin*  s  m.aximum  principle 
then  allows  one  to  solve  the  S.O.C.  problem  by  solving  a 
two-point  boundary-value  problem. 

3.   Gradient  Technique 

This  technique  is  a  successive  approximation  pro- 
cedure in  which  an  iteration  is  made  on  the  control  vari- 
able u(t)  so  as  to  improve  the  function  to  be  minimized 
at  each  new  iterate,  subject  to  the  constraining  differ- 
ential equations. 

Consider,  for  exam.ple,  those  problems  in  which  it 
is  required  to  extremalize  performance  indices  of  the 
form 
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J(x,u,t)  =      f      (t>    (x.u.t)  dt  (2.15) 

•'o 


subject  to  the  constraint  in  the  form  of  the  first  order 
differential  equations 

0i  =  Xi  -  "^^(x.H.t)  =  o  •    (2.16) 

It  is  required  to  find  u*(t)  which  at  the  same  time  ex- 
tremalizes  J  and  satisfies  equation  (2.16).  The  overall 
algorithm  which  simultaneously  satisfies  both  conditions 
gives  (i+l)st  iteration  of  u  in  terms  of  the  ith  iteration: 

m 


'-^  =  "^tK[^.   .,  l-^] 


(2.17) 


where  m  denotes  the  number  of  equations  9   and  K  and  A  are 
chosen  arbitrarily.   It  was  seen  in  the  section 
that  the  maximum  principle  applied  to  the  above  set  of 
equations,  resulted  in  a  coupling  equation,  coupling  the 
state  vector  and  the  adjoint  vector  via  the  control  vari- 
able u.   This  result  therefore  can  modify  equation  (2.17) 
to  yield  the  following  algorithm  which  approximates  some  of 
the  characteristics  of  (2.17), 

pi+1  _  pi  +  K  [  pi  (tf)  -  P^(tf)]  (2.18) 

This  algorithm  was  used  in  the  problem  solved  in  section 
B.l.  This  is  an  excellent  approximation  to  (2.17)  which 
is  most  useful  on  a  repetitive  analog  computer. 
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4  .   Quasllinearization  Technique 

This  technique, developed  by  Kalaba,  Bellman  and  others 
is  particularly  suited  for  the  type  of  optimization  problems 
being  considered  here.   It  provides  an  effective  computa- 
tional tool  for  a  wide  class  of  non-linear  two-point  and 
multi-point  boundary  value  problems  (15.18.19).   The 
technique  consists  of  producing  a  sequence  of  approxi- 
mating functions  xO(t).  x^^ ^ (t ) . . . , x^") (t )  which  converge 
rapidly  to  the  solution  of  the  original  problem.   In  (14.20) 
it  is  shown  that  this  convergence  is  quadratic,  i.e., 
asymptotically  the  number  of  correct  digits  in  each  approx- 
imation is  doubled. 

To  illustrate  the  technique,  consider  a  2N-dimensional 
dynamical  vector  equation: 

^=     |(>^)  (2.19) 

subject  to  the  boundary  conditions 

Xi(o)  =  bi      i  =  1.2, ..,,N 

(^    \        u  (2.20) 

Xi(tf )  =  bj_     i  =  N+1,...  ,2N 

Assuming  that  the  nth  approximation  is  known,  find 
the  (n+1)^  approximation  as  the  solution  of  the  linear 
two-point  boundary  value  problem.  To  do  this  take  a   differ- 
ential of  (2,19)  and  obtain 
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A7-    ^'(^) 

.   Ax                          (2.21) 

bx 

or   x^k+l  -  x^l^  = 

^  ^.ilill^   (   (k.l).,Jk))  ^ 
^'                 '-                  (2.22) 

i  =  1.2 

From  (2.19)  above 

ii^=  ^l(xK) 

(2.23) 

Therefore 

C  k+1  _  ,  /^ks 

v^    ^n(x'<)  ,   n,4.i  \    ,.  .  . 

J.,       a i.  •       ■  ^^  .      'Tn     ~x^^^'  v^.24) 

m=l     °  ^m  m 

Mote  that  the  terms  on  the  right-hand  side  of  equation 
(2,24)  are  the  first  two  terms  in  the  power  series  expan- 
sion of  the  function  ^i(x    ')about  the  point  x*^. 

The  computational  scheme  consists  in  the  numerical 
production  of  a  particular  solution  and  two  independent 
homogeneous  solutions  of  equation  (2.24)  on  the  interval 
(o.tf)  and  the  determination  of  the  constant  multiplier 
of  the  homogeneous  solutions.   The  solution  may  be  repre- 
sented in  the  following  form: 

x^-'Ut)   =  p  (t)  +    E  Ci  hi(t)  (2.25) 

i=l 
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where  the  vector  function  p(t)  is  the  particular  solution 
of  equation  (2.24),  i.e.: 

P=  ^i^^*")  ^  ^     ^HJ  (P-Xm^^))       (2.26) 

Subject  to  the  initial  conditions 

p{o)  =  0 

and  the  vector  function  hj^(t)  is  the  solution  of  the  homo- 
geneous form  of  equation  (2.24),  i.e.: 

m-1      o^m 

The  initial  vector  has  the  ith  row  unity  and  the  others 
zero.   These  vectors,  p(t)  and  hi(t),  are  all  considered  to 
be  determined  on  the  interval  o  £  t  <  t£ . 

The  solution  of  equation  (2.24)  i.e.  equation  (2.25) 
together  with  the  boundary  conditions  in  <  quation  (2.20) 
leads  to  the  following  linear  algebraic  equations: 

x't^-^(o)  =  b^=  p(o)  +   E  cjhj(o)  ,   -L=i,^,   -   (2.28) 
*       j=l 

and   xl^"^l(tf)  =  b,=  pjtf)  +   E  c^h^^itf),  ^'-^^\y-,^^    (2.29) 

In  equation  (2.28)  since  p(o)  and  hi(o)  are  known 
the  scalars  c^  appearing  in  equation  (2.28)  and  (2.29) 
are  determined  in  terms  of  the  boundary  conditions. 
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Similarly  the  function  x    is  determined  from  the 

— k+1 
known  values  of  the  function  x   .   In  this  manner  a 

large  system  of  non-linear  differential  equations  subject 

to  the  boundary  values  may  be  resolved. 

B.   Formulation  of  Specific  Problems 


1 .   Cost  Function  1 


J  =  1/2  J         Pc^(t)dt   subject  to  the  constraints  that 
o  n(tf)  =  ano 

TPBVP:   By  using  the  statement  of  the  maximum  principle 
as  in  A. 2,  write  the  Hamiltonian  as: 

H(n,c,pi.P2.  p  ,t)   -   1/2  p^^   +  P2(t)n+P2(t)c 

■=   1/2  Pc^  +  Pi(t)[  ;\(  P^-  B)n^   /.c] 
+  P2(t)[  ;f  n  -  Ac]        (2.30) 

Then  by  (2.8) 

aH(n,c,p, ,pp,  p  ,t) 

P,  (t)  = -^— ^ (2.31) 

-I-  dn 

and 

dH(n.c,Pi,P2,  P  ,t) 
P2(t)  = ^-^ (2.32) 

From  equations  (2.31)  and  (2.32)  get  the  auxiliary  equa- 
tions 

Pl(t)  =  -  ^  (p-/3)Pi(t)  -  4  P2(t)  (2.33) 
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and 

P2(t)  =  -  /.P^d)  +  ;.P2(t)  (2.34) 

Setting  -^  =  0,  yields 

-^c  =  !l!  (2.35) 

A 

The  following  system  of  equations  are  therefore  obtained. 
r^  =  -i   Pc  -   B  )    r\  -^     /.c  (2.36) 

6  =  -^  n  -  AC  (2.37) 

Pi  =  "  A  ^  ^c  -  ^  )  Pi  -  T  ^2  (2.38) 

P2  =  -  aP^  +  /.P2  (2.39) 

P,n 
Pc  =-  A  (2.40) 

with  boundary  conditions: 

n(o)  =  no 
C(o)  =  c^ 
n(t^)  =  an^  (2.41) 

P2(tf)  =c.C 


ly 


Simulation  of  problem  1 

The  simulation  of  the  problem  as  described  above,  for 
nuclear  reactor  system  dynamics,  was  performed  with  an 
Applied  Dynamics  analog  computer.   A  detailed  computer 
wiring  diagram  used  for  programming  the  computer  is 
given  in  Figure   41  .   Table  1  gives  a  list  of  system 
constants  and  the  initial  conditions  used.   The  system 
equations  described  in  Section  1,  including  the  constants 
of  Table  1  are  given  below. 

Problem  1; 

n  =  1000  p^n-6.4n  +   ,   Ic  (2.42) 

c  =  6.4n  -  .Ic  (2.43) 

F^  =-1000  ^cPi  +  6.4  P^  -  6.4P2  (2.44) 

P2  =  .1  Pi  +  .1  P2  (2.45) 

Pc  =-1000  nP^  (2.46) 

With  boundary  conditions: 
n(o)  =  5.0 
c(o)  =  64n^ 
n(tf )  =  5.0a 
P2(tf)  =0,0  (2.47) 
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TABLE  I 


N'JCLEAR  REACTOR  CONSTANTS  FOR  A  LUMPED  SYSTEM  AND  INITIAL 
CONDITIONS 


Nuclear  Reactor  Constants 


Decay  constant,  /.  OclO  sec~^ 

Average  neutron  life  time, A  10~'^  sees 

Delayed  neutron  fraction,^  0.0064 


Initial  Conditions 


•^'eutron  density  or  power  level,  n  5K.W 

Precursor  population,  c  64n_, 
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Simulation  of  the  system  of  equations  on  an  analog 
computer  requires  amplitude  scaling.   This  scaling  is 
necessary  because  all  voltages  must  be  kept  within  +  100 
volts  in  order  to  prevent  overloading  of  the  machine 
components.   This  overloading  will  lead  to  inaccurate 
results.   The  basic  set  of  scaled  equations  are  rewritten 
here  from,  Appendix  A  for  convenience.   The  bars  indicate 
scaled  quantities. 

Problem  1 : 

N  =-100  N^PjL  -  6.4  N  +  10  C 

C  r  0,064N  -  0.1  C 

P  =   100  Np2  +  6.4Pj^-6.4P2 

P2  =  -  0.1  P^  +  0.1  P2 

Pc  =-0.1  N  P^ 
With  boundary  conditions: 


K  =  ^' 

,05 

Co  =  0. 

,032 

N(tf)   = 

=  aNo 

P2(tf) 

=  0. 

0 
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Typical  solution  procedure:   With  the  equations  scaled,  the 
modified  gradient  rricthod  will  be  used  for  optiaization  on 
the  analog  computer  and  the  results  verified  with  the  num- 
erical technique  of  quasilinearization. 

The  gradient  technique  is  preferred  as  opposed  to 
the  tria]  and  error  procedure  which  can  be  quite  cumber- 
some and  extremely  tim.e  consuming,  particularly  for  high 
order  systems,  A  brief  description  of  this  technique  is 
given  in  Section  A. 4. 

The  equations  used  for  simulation  are  on  page  21, 
and  the  circuit  diagram  is  shown  in  Figure  ^1 , in  Appendix 
A.   The  neutron  density  or  the  power  level  is  forced  to 
different  desired  values  ranging  from  n^   =  2n  to  n^   =  lOno 
in  one  second.   Figure  1  shows  the  neutron  density  plots 
and  the  reactivity  function  plots.   Table  2  shows  the 
initial  values  of  the  adjoint  variables  and  the  cost 
function  to  attain  various  peak  power  levels.   It  is  ob- 
served that  n(tf)  ^0  for  all  t  >  tf.   To  accomplish  this, 
the  reactivity  will  have  to  be  programmed  differentially 
for  all  t  >  tr.      This  can  be  easily  determined  from  the 
kinetic  equations  (  10,11) . 

The  cost  function  is  plotted  versus  the  peak  power 
attained.   This  is  shown  in  Figure  2.   The  graph  indicates 
increase  in  cost  with  increase  in  power  level  and  this  in 
general  is  to  be  expected. 


I. 


,4         .6 

TIME  (Sec.) 


Figure  1.   Power  level  N  and  reactivity  Q^  versus  tirr.e 
for  cptirral  control  problori',  1.   (Nuclear 
Reactor  System) 


24 


in 
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DESIRED  POWER  LEVEL,  N^j 
Figure  2.   Cost  function  versus  desired  pov;e: 


level 
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TABLE  2 


INITIAL  CCMDITir^'S  PN  THE  ADJCI^^T  VARIABI  ES  AND  THE  COST 
J.  TO  ACHIE\'E  DESIRED  PCWER  LF\'FI,S 


N^ 

J 

PlO 

P20 

2no 

0.454x10'"^ 
1.263x10'"^ 
1,51    xIO"^ 
1.945x10"^ 

0.629x10"^ 
0.631x10"^ 
0.552x10"'^ 
0.457x10"^ 

0.408x10"^ 
0.268x10"^ 
0.235x10"^ 
0.176xl0"^ 

The  gradient  technique  is  studied  in  greater  detail 
for  only  one  case,  namely  that  of  n^  -   lOn^^,   A  procedure 
is  developed  whereby  changing  P^ ,  a  search  for  P^*  which 
produces  the  desired  optirral  control  u*  (reactivity  in 
this  case)  and  trajectory  n*  can  be  made.   Actually,  if 
all  P^  values  are  scanned,  much  more  general  information 
concerning  the  control  of  the  system  may  be  obtained. 
In  particular  the  reachable  set  of  the  system  may  be 
determined.   The  reachable  set  RCt.n^)  is  the  set  of  all 
states  which  may  be  obtained  at  time  t  (starting  from  n^ 
at  time  t  =  0)  by  means  of  admissible  controls.   Thus 
R(t,nQ)  determines  the  extent  of  possible  trajectory 
with  admissible  controls.   By  taking  a  suitable  set  of 
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Pi's  enough  boundary  points  may  be  obtained  to  define 
R(t^,nQ),  (tf  is  a  particular  value  of  t,  one  second  in 
this  case).   This  can  easily  be  seen  in  Table  3.   Run 
C-12A  for  instance  is  plotted  in  Figure  (3).   Converg- 
ence of  the  trajectory  to  the  desired  one  is  shown. 
Only  the  last  four  iterations  are  plotted.   Correspond- 
ing control  functions  are  also  shown. 

Updating  Circuit :   The  equation  for  the  iterative  pro- 
cedure used  are  obtained  from  Equation  (2.18),  and  the 
computer  circuit  is  given  in  Figure  44   ^  Appendix  A. 
The  equations  are: 

P^J^^  =   PJ(0)  1  Kii[N^(t  )  -  Nd(tf)] 

-  ^^4^^^f^-P2d(^f)J  ''-''' 

p^J-^l  =  P2^(0)  t   K2i[  ^'^(t)  -  Nd(tf)] 

1  K22[  P^Ctf)  -P2d(tf)-]  (2.49) 

The  above  search  algorithm  seem.s  to  be  the  logical 
choice  since  it  is  desired  to  satisfy  the  following  end 
conditions  : 

N(tf)  =  lOn^ 
P2(tx)  =  0 
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The  last  updated  values  of  the  P^^'s  for  which  the 
above  conditions  are  satisfied  are  held  on  the  n  analog 
memories  (number  two  integrators  of  the  track-and-hold 
system.  Figure   44  ).   The  size  K^x-  Kx2.  K2X ,  K22  of 
the  increments  8 P^  are  selected  by  the  potentiometers. 
The  convergence  of  the  iterations  is  very  sensitive  to 
the  values  of  K's.   Larger  values  of  K  result  in  large 
overshoots  in  the  desired  trajectory  (see  comment  column 
in  Tablel2).   Smaller  values  result  in  an  increased 
number  of  iterations  to  converge  -  Figure  4  .   The 
tolerance  on  the  end  conditions  during  these  tests  are 
of  the  order  of  a  few  volts.   No  attempts  are  made  to 
obtain  high  accuracy  in  the  terminal  conditions  as  this 
would  require  use  of  more  analog  and  logic  elements  and 
this  is  not  warranted  at  this  stage.   Figure  3   shows 
the  sequence  of  trajectories  converging  to  the  optimal 
one.   In  a  typical  case  with  P^(o)  far  from  optimal, 
convergence  was  obtained  in  twenty-nine  trials.   Total 
time  for  convergence  was  60  seconds.   A  graph  of  number 
of  iterations  for  convergence  versus  Piq^  is  shown  in 
Figure   4  for  various  step  sizes.   Note  that  when  P.„^ 
is  very  close  to  the  optimal  value,  the  number  of  itera- 
tions for  convergence  is  infinite;   Figure   5 
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shows  the  plot  of  the  maximum  value  of  n(neutron  density) 
and  the  value  of  Poq  versus  the  iteration  number  in  a 
particular  iteration  cycle.   With  successive  iterations 
the  value  of  P20  decreases  until  an  optimal  value  is 
obtained,  i.e.,  value  for  which  n(tf)  -   ICn^  and  P2(t£) 
=  0.   Note  that  in  C-12A  since  one  is  closer  to  the 
optimal  case,  the  change  in  P20  is  gradual  as  also  with 
maximum  n.   Moving  away  from  optimal  setting  of  Poq*  ^^® 
change  in  P20  is  much  more  and  this  is  similarly  trans- 
lated to  maximum  n. 

The  steps  in  the  iterative  procedure  are  sufficiently 
complex  that  it  is  wise  to  break  the  programming  down  into 
a  series  of  sub-problems. 

An  iterative  procedure  consists  of  the  following  steps 

(a)  Choose  an  initial  Pj^J(o)  to  be  the  same  as 
that  obtained  by  trial  and  error.   Let  P2'^(0) 
take  on  arbitrary  values. 

(b)  Motice  response  of  the  power  level  (n).   This 
will  help  to  determine  sign  of  8P^.   During 
the  iteration  process,  the  addition  is  per- 
formed by  the  analog  m.emoiry  so  that  the  new 
value 

^2    "■  20       2 
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where 
SP2^  =  t  K21  rnJ(t)  -  n^(tf)  1  K22rP^(tf)  "  P2d(tf 

can  be  transferred  to  the  P20  memory,  when  up- 
dating of  this  parameter  is  required.   The 
values  of  K22  and  K22  are  chosen  to  satisfy  the 
tolerance  and  oscillation  conditions  mentioned 
earlier. 

(c)  Now  reverse  the  above  procedure.   Choose  an 
initial  P2-^(0)  as  obtained  by  trial  and  error. 
Let  P,-^(0)  take  on  arbitrary  values. 

(d)  Observe  response  of  power  level  to  determine 
sign  of  ^  Pj^ .   During  the  iteration  process 
addition  is  performed  according  to  the  follow- 
ing equation 


J  +  1  - 


10 


+  ^P, 


Procedure  in  (b)  is  repeated  for  suitable 


selection  of  K^^^  ^^^   ^ 


12' 


(e)   Having  determined  the  signs  of  hP^    and  i^P^, 
P2^-^{o)  and  P2"'(o)  are  now  assigned  arbitrary 
values  and  convergence  to  the  true  values  is 
obtained  by  observing  the  neutron  density  re- 
sponse during  each  operation  cycle  to  the  end 
of  the  iterations.   The  computer  circut  dia- 
gram for  updating  P-"  is  shown  in  Figure 
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Figure  3.   Power  level  and  associated  reactivity 
time  plots,  leading  to  the  desired 
trajectory  using  gradient  technique. 
(Muclear  Reactor  System) 
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Figure  4.   Number  of  iterations  for  convergence  versus 
PlO 
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Problem  1  ^la  Quasilinearization  Technique:   The  quasi- 
linearization  technique  is  described  on  page  14  and  will 
now  be  applied  to  this  problerr. 

For  convenience  rewrite  equations  (2.42-2.45)  and 
modified  by  (2.46)  as : 

(2.50) 


n  =  5xl0^pin2-6.4n+0.1c 

c  =  6.4n  -  0.1c 

P^  =  5xl0^Pi2n+6.4P^-6.4P2 

?2   =-0aPi+0.1P2 


(2.51) 
(2.52) 
(2.53) 


Linearizing  the  above  equations,  as  on  page  15,  there  results 


pk+l  =  nk  +  5x10^ 


'pik2nk(n^^l-nk).(n^)2(p^k.l.p^k)j 


6.4(n'^^i  =  nk)  +  0.1  {c^^^-c^) 


But 


n*'  =  5xl0^pi''(n'')^  -  6.4n''  +  0.1  c^ 


Therefore 


k+1 


=  5xl0^pi'^(n'^)^-6.4n'^  +  O.lc"^  +  5x10' 


k^  k 
Pi  2n 


I    k+l   kv  ,  / 
(n   -n  )  +  ( 


n'^)'(Pl'^^-p/)]  -  6.4(n^*l-n' 


+  0.1  (c'^^l-c'') 
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Hence 


-  6.4n''^^  +  O.lc'''^^  (2.54) 


Mow  let  n  =  x^,  c  =  X2,  p^  =  X3  and  p2  =  x^,  and  let 
k  =  o  and  this  leads  to 

xi^   =  10^xOx^l-106x3O(xiO)2+5xl05x3^(x^°)^ 

-  6.4x^^+0. 1X2^  (2.55) 

Linearizing  Equation  (2.51)  there  results 

ek+1  ^  -k  +  6.4(n'^^l.n'^)  -  OA{c^^^-c^) 

=   6.4n'<-C.lc»^+6.4(n''"'^-n'^)  -  0.  Kc'^^'^-c'^) 

or   c'<^i  =  6.4nk^l  -  Clc^+l  (2.56) 

let  k  =  o,n  =  XI,    etc,  and  get 

X2^  =  6.4x^1  -  0.1x2^ 

Repeating  above  procedure  for  equations  (2.52)  and  (2.53) 
obtain : 

^3!  =  10^xi°(x30)2  .  106x^°X3Ox3l»5xl0^(x3C)2x/ 

+  6.4x3!  -  6.4x4!  (2.57) 
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and  X 


•  .1 


4-^  =  -O.lxol  +  0.1: 


(2.58) 


The  above  equations  tray  be  written  in  the  following  form: 


^1^   =  b^  +  aiix/  +  ai2X2^  +  ai3^3^ 


i  + 


<3^  =  b3  +  331X^1  +  333x3^  +  334x4! 
<4'^  =  343X3^  +  344X4! 


In  vector  notation  x 
where 


•N-1   =A(t)    x^^^l.   b(t)    ^^1 


(2.59) 

(2.60) 

(2.61) 
(2.62) 

(2.63) 


X   = 


b  = 


and   the   components    of   the  A  inatrix   and   b  vector    are: 


106(xi°)2   X  3° 


b2     =  b4  =  0 


b3     r   10^i°(x3°)2 


a^i   =   10^   X   1^X3-6.4 


^12  =  0.1 

ai3  =  5xlO^(xi°)^ 
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^14  -  ^23  =  ^24  =  332  =  ^41  =  844  =  0 

321  =  6.  4 

322  =0.1 

331  =  -5xl0^(x3°)^ 

333  =  -10^x1^x3°  +  6.4 

334  =  -6.4 

343  =  -0.1 

344  =0.1 

It  should  be  noted  here  thst  in  vector  equ3tion  (2.63), 
the  coefficients  3re  time  vsrying.   Equ3tion  (2.63)  m3y 
also  be  rewritten  3S 

x^  =   L  a/(t)xJ  +  b^(t).  i  =  l,2.._,n.     (2.64} 

If  p(t)  is  the  psrticulsr  solution  of  equation  (2.63).  then 
P^  =   E    3/(t)pJ  +  bi(t)  (2.65) 

subject  to  the  initial  conditions  p(o)  =  0 
or 

p^  =  h^  +  a^^n  +  ai2P2  ^   ^13^3  (2.66) 

P2  =  ^21^1  ^  a22P2  (2.67) 


38 


P3  =  ^"3  -^  331  ""  ^33  ^^-^^^ 

P4  =  a43P3  ^  a44P4  (2.69) 

Appendix  B  shows  incorporation  of  equation  (2,65)  in  the 
main  programme. 

The  homogeneous  solution  of  (2.64)  is  the  solution  of 
equation 

i.  =   £  a.^(t)  xJ  .  i  =  1.2... .,n  (2.70) 

^    j=l   ^ 

or  3c  =  A(t)x  (2.71) 

Pontryagin  (20)  points  out  that  the  solution  of  the  above 
equation  with  time  varying  coefficients  is 

x(t)  =   0  (t  tjx(to)  (2o72) 

where 

i(t)  =  A(t)  :'(t),  -(o)  =  I  (2,73) 

where  c(t)  is  an  nxn  matrix 
A(t)  is  an  nxn  matrix 
I  is  a  unit  nxn  matrix 
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The  complete  solution  is  therefore 

x(t)  =  4)  (t)x(o)  +  p  (2.74) 

The  following  initial  and  final  conditions  were  used 

n^  =  5  k.w  or  0.5  n  (t=tf)  =  lOn^  =  5 

32  P2(t=tf)  =  0 


10 


0.5  X  10-^  Pi(tf )  =  0,30x10=^    ^^-"^^^ 


P20 


=  0.18  X  10= 


The  initial  guesses   for  the  ad  joint  variables ,  i.e.  p^Q 
and  P20  were  taken  from  the  analog  computer  results 
described  earlier.   The  above  case  is  for:   neutron  life 
time  A  =  10;^^3  . 

The  numerical  results  of  the  first  control  problem 
may  be  summarized  in  the  graphs.   The  graphs  also  indicate 
comparison  with  the  analog  computer  results,  which  in 
general  is  very  good.   Two  cases  were  run  -  average 
neutron  life  time,  y\  =  10~^secs  and  A  =  10~^secs.   A 
typical  table,  i.e.,  Tables,  shows  the  results  of  the 
second  case  only.   Comparison  with  analog  computer  results 
is  the  subject  of  Figure  (6). 
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The  initial  and  final  conditions  used  in  Figure  7 
are: 

Pj^(t£)  =  0„0000020,  P2(tf)  =  0.0 

n(o)   =  0.5,-  c(o)  =  32  and  a  =  10"^ 

Figure  8  shows  power  level  versus  time  plot  of  run 
C-12A  of  Table  12.  The  initial  conditions  on  the  adjoint 

variables  in  this  run  are:  p^Q   =  0^54x10   and  P20  = 

—5 
0<,18xlO   .   Convergence  in  this  case  is  obtained  in  the 

second  trial. 
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Corrparison  of  analog  computer  re- 
sults with  those  obtained  by  quasi- 
linearization  (Nuclear  Reactor 
System) 
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Power  level  versus  time  (  A 
Comparison  of  analog  computer  results 
with  those  obtained  by  quasilinearization 
(Nuclear  Reactor  System) 
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2.      Cost   Function  2 


=  1/2     f^[Pc2  +    a(n-n^)^    1 


dt 


The  two-point  boundary  value  problem  resulting  by 
the  application  of  the  maximum  principle  is: 

(2.76) 

c  =  6.4n  -  0.1  c  (2.77) 

P^  =   a(n-nd)  -  1000  p^P^   +  6.4  p^-6.4P2     (2.78) 
P2  =-0.1  P^   +  0.1  ?2  (2.79) 

Pc  =-1000  nPi  (2.80) 

with  boundary  conditions: 

n(o)  =  5-0 

c(o)  =  64n^  (2.81) 

n(tf )  =  5.0a 

P2(tf)  -  0.0 

The  cost  function  of  this  problem  may  be  written  as 
the  sum  of  the  two  cost  functions  as  shown. 

J  =  1/2   j  [  ^c2  +  a  (n-nj)^j  dt   =  Ji  +  ^2 
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where  J^  =  1/2   J  P  c^  (it 
o 


and  may  be  thought  of  as  the  cost  of  the  control  action  and 


^f         9 
J2  =1/2  j     a   (n-nd)^dt 
o 


is  the  cost  of  the  deviation  of  the  output  of  the  dynamical 
system  from  the  desired  value. 

Simulation  of  Problem  2 

This  problem  reduces  to  problem  1  for  a  =  0,   A  graph 
of  the  cost  function  for  various  values  of  a  is  shown  in 
Figure   9  ,   Neutron  density  and  reactivity  function  plots 
are  shewn  in  Figures  10,11  for  ae  (0,20)  with  the  final 
time  tf  fixed  at  1  sec.   It  is  clear  from  these  figures 
that  for  increasing  values  of  a  ,  the  rate  of  control  rod 
withdrawal  is  also  increased.   However,  in  order  to  limit 
the  power  level  to  a  certain  desired  value  the 
control  rod  is  inserted  back  into  the  reactor  at  a  faster 
rate,  thus  shifting  the  reactivity  function  peak  to  the 
left. 


47 


Table  5  shows  the  peak  power  attained,  a   and  the  cost 
functions  Jj^,  J2  and  J.    With  appropriate  scaling  J2 
is  given  by: 


10"3  _  tf      o 

-T—  a  j   0.1  (N-Nd)2dt 


J2  =  5  X  lO""^   J2 
where 


J2"  =  a  J   O.l(N-Njj)  dt 


From  the  graph  of  Ji  versus  O    (Figure  12)  it  appears 
that  the  cost  function  increases  with  increasing  a   .      This 
appears  reasonable  as  increased  a  would  decrease  the  error 
between  n  and  n^   and  hence  more  control  would  be  needed. 
The  total  cost  function  increases  rapidly  with  an  increase 
in  a  .  as  is  apparent  from  the  graph  of  the  total  cost  J 
versus  N,  -  the  peak  power  attained  (Figure  13).  and  also 
from  graph  of  J  versus  a  ,   This  result  is  again  quite 
logical  as  one  would  expect  to  pay  a  higher  cost  in  crder 
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Figure  9, 
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i"otal  cost  function  versus  peak  power 
desired  (control  problem  2) 
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TABLE  4 

COST  FUNCTIONS  AS  A  FUNCTION  OF  WEIGHTING  0  AND  POWER 
LE>/EL  N^  FOR  PROBLEM  2 


2Ji 


a  =z  5.0 


J=Jl+J2 


2no  0,456x10"^ 
5no  0.696x10"^ 
6no  1.323x10"^ 
lOrio  1.914x10"^ 


-5 


0.0149x10 
0.527x10"^ 
0.823x10"^ 
2.88x10"^ 
a  =  10 


0.243x10"^ 
0.875x10"^ 
1.484x10""^ 


3.84x10" 


,-5 


2nQ  0.456x10  ^ 

5no  0.696x10"^ 
.-5 


0.0297x10 


,-5 


0.258x10"^ 


IOpq  2.058x10"^ 


4.95x10  ^ 
a  =  20 


6,0x10"^ 


2no 

0,456x10"^ 

0.0595x10 

^% 

0.745x10"^ 

1.92x10"^ 

6no 

1.445x10"^ 

2.81x10"^ 

IOHq  2.126x10 


-5 


9.2x10 


0.288x10 
2.30x10"* 
3.53x10"' 


10.26x10 
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Figure  12.  Cost  function  ? 
movement  versur> 
pi- obi  em  2) 


iSociated  with  contio].  rod 


M.- 


10  n.,  (control 
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Figure  13.   Total  cost  function  versus   ^.  ,  for  control 
problerr.  2  end  n^   =:  lOn^^ 
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3.   Cost  Function  3 


f  _  9  .  ^  / ^2 


J  =  1/2   J   Pc^  +  a  (n-n^) 
o 


dt 


p  Q  =  a  +   A(n-nci) 

The  two-point  boundary  value  problem  resulting  by 
the  application  of  the  maximurr.  principle  is: 

n  =  1000  an  +  1000  A(n-n^)  -  6.4n  +  0.1c       (2.82) 
c   =  6.4n  -  0.1c  (2.83) 

(2.84) 


A  =  0.  a  =  0 


P^  =  -  (a2  +  a  )(n-n^)  -  aA  -  1000  P^ (a+2An-Anj) 


+  6.4?^  -  6.4P2 


P2  =  -  O.lPi  +  0.1  P2 


(2.85) 
(2.86) 


P^  =  -  A(n-nj)^  -  1000  np^(n-n^)  -  a(n-nd)     (2.87) 


P4  =  -a  -  A(n-nd)  -  lOOOnpj^ 


(2.88) 


Pc  -  a  +  A(n-nd)  (2.89) 

The  boundary  conditions  are: 

n(o)  =  5.0  P2(tf)  =  0 

c(o)  =  64nQ  n(tf)  =  5.0  a             (2.90) 

P3(o)  =  P4(o)  z:  0  P3(tf)  =  P4(tf)  =  0 
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Simulation  of  Problem  3 

The  equations  used  for  simulation  are  given  in 
Appendix  A,  and  the  procedure  used  is  similar  to  that 
followed  in  Problem  1.   The  circuit  diagram  is  shown  in 
Appendix  A.   The  trial  and  error  procedure  is  adopted  in 
the  sense  that  the  constants  "a"  and  "A"  are  varied  to 
obtain  the  desired  power  level  response.   Figure  14  shows 
the  variation  of  "A"  with  "a"  to  achieve  different  power 
levels  ranging  from  Sn^  to  Sn^.   Corresponding  to  these 
power  levels,  the  desired  control  is  obtained.   The  initial 
conditions  of  P^  and  P2  and  the  weighting  factor  a    are  now 
varied  to  satisfy  the  end  conditions  on  the  adjoint  vari- 
ables, namely  P2(tf)  =  Paltf)  =  P4(tf)  =  0. 

It  is  observed  in  this  investigation  and  justified  by 
solution  of  equation:   N  =  AN(N-N(j)  -  0.64N  +  C  that  the 
desired  power  level  is  reached  only  for  an  infinitely 
large  value  of  A,  Figure  (15).   The  power  level  response 
also  suggests  a  strong  dependence  on  the  desired  level. 
These  observations  therefore  suggest  that  the  form  chosen 
for  this  control  will  not  yield  a  desirable  closed  loop 
control.   The  form  that  accomplishes  this  is  given  in 
Problem  4. 


Figure  14.   V/ariation  of  parameters  in  the  con- 
trol law  for  s.o.c.  problem  (i.e., 
problem  3)  and  for  various  values 

of  V.^ 
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8n^ 


26n, 


4n, 


VALUE  OF 


Figure  15.   Peak  power  level  reached  versus  A  (control 
problem  3) 
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4„      Cost   Function  ^ 


J  =  1/2         l^    [Pc^  ''    ^  ("-"d)^J  dt 

t£ 

p^  =  A(n-n^)    +   B       J  (n-nd)dt 

o 

Let  X  =   J  (n-n^^)dt  Hence  x  =  (n-n^).  and  p^  =   A(n-nc|)  +  Bx 

o 


The  two-point  boundary  value  problem  resulting  by  the  appli- 
cation of  the  maximum  principle  is: 

n  =  1000  A  n(n-nd)  +  1000  Bxn  -  6.4n  +0.1c      (2.91) 

c  =  6.4n  -  0.1c  (2.92) 


X  =  n-n^j 


(2.93) 


A  =  0  .  B  =  0  (2,94) 

P2  =  (A^+a  )(n-n(j)  -  ABx  -  1000  P^  F  2An=Anj+BxJ 

+  6.4?^  -  6.4P2  -  P3  (2.95) 

P2  =  -O.IP^  _  O.IP2  (2.96) 

P3  =  -B^x  -  1000  BnPi  -  AB(n-nj)  (2.97) 

P4  =  -A(n-n^)^  -  lOOOnP^(n-n^)  -  Bx(n-n^)  (2.98) 

P5  =-Bx2  -  1000  xnP^  -  A  x  (n-n^)  (2,99) 

p^   =  A(n=nd)  +  Bx  (2.100) 


59 


(2.101) 


The  boundary  conditions  are: 

n(o)  =  no 

c(o)  =  Co 

x(o)  =  0 

n<^(tf )  =  n^d  =  ano 

P2(t£)  =  0 

P3(tf)  =  0 

P4(tf)  =  0 

Psdf)  =  0 

P4(o)  =  0 

Pp^(o)  =  0 

Simulation  of  Probletr  4 

It  was  seen  in  Problem  3  that  the  power  level  re- 
sponse was  highly  dependent  on  the  setting  of  the  desired 
value  of  the  neutron  density.   In  order  to  remove  this 
restriction  a  control  function  of  the  type  shown  above, 
namely  a  combination  of  proportional  and  integral  control, 
was  tried.   A  simple  transfer  function  analysis  indicates 
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that  such  a  system  would  be  stable  and  the  response  would 
be  dependent  on  the  coefficients  A  and  B. 

The  equations  used  for  simulation  are  given  on  page 
63  and  the  procedure  used  is  similar  to  that  followed  in 
Problem  3.   Circuit  diagram  is  shown  in  Appendix  A.   Figure 

16  shows  the  power  level  and  the  reactivity  function 
response  for  values  of  n^j  ranging  fron  2nQ  to  lOn^  with 
the  final  time  fixed  at  one  sec.   Note  in  this  case  that 
the  power  level  rise  is  much  faster,  even  though  it  settles 
down  to  the  desired  value  in  fixed  time  tf .      Corresponding 
rod  withdrawal  is  also  much  faster  in  the  initial  stages 
and  then  gradual  before  it  reaches  a  constant  value,  to 
maintain  the  above  level.   Figure  17   shows  variation 
of  cost  function  versus  peak  power  attained.   A  graph  is 
also  plotted  showing  behavior  of  constants  A  and  B  with 
respect  to  peak  power  attained  (Figure  18).   It  is  clear 
from  this  graph  that  there  is  no  value  of  optirrum  A  and 
B  which  is  independent  of  the  desired  power  level.   This 
would  imply  that  the  control  parameters  A  and  B  should 
not  be  fixed  but  should  be  varied  as  a  function  of  the  de- 
sired final  state.   Hence  one  would  conclude  that  some  sort 
of  adaptive  control  may  be  required.   A  control  system  is 
said  to  be  adaptive  if  it  is  capable  of  changing  its  con- 
trol object  in  such  a  way  as  to  operate  at  all  times  in 
an  optimal  or  nearly  optim.al  fashion. 
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Figures  19   and  20  indicate  the  time  response  of 
the  auxiliary  functions.   Since  these  functions  are  highly 
unstable,  they  are  very  sensitive  to  amplifier  drift. 
Because  of  this  reason  it  is  very  difficult  to  satisfy  and 
maintain  the  end  conditions,  namely,  P2(t£)  =  Po(tf  )  = 
P4(tf )  =  P5(t£)  =:  0,  and  n(tf)  =  an^. 


TABLE  5 


COST  FUNCTION  AS  A  FUNCTION  OF  POWER  LEVEL 
N^  AND  CONSTAhTTS  A  AND  B  FOR  PROBLEM  4 


(NUCLEAR  REACTOR  SYSTEM) 


^1 


lOno 

2.123 

42.12 

3.178 

1.2C1  X  10"^ 

6^0 

2.709 

55.50 

5.667 

0.495  X  10~^ 

5no 

3.695 

79.78 

11.34 

0.294  X  lO"^ 

2no 

3.836 

128.4 

" 

0.105  X  10-5 
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Nd=IOno 


A  .6 

TIME  (Sec.) 

Figure  16.   Pov/er  level  ^J  and  reactivity  c^  versus 
time  for  specific  optimal  control  pro- 
blem 4  (Nuclear  Reactor  System) 
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Figure  17.   Cost  function  versus  peak  power 
attained 
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Figure  18.   Nj  versus  ^   and  B 
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Figure  19.   Auxiliary  functions,  P2 ,  P3,  F4  and 
PC)  versus  time  for  control  problem 
4;  N(j  -  lOn^j  and  tn^    (Muclear  Reactor 
Systen-i) 
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Figure  20.   Auxiliary  functions  P2,  P3 ,  P4  and  P^ 
versus  time  for  control  problem.  4; 
Nd  =  5n^^  and  2no  (Nuclear  Reactor  System) 
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Discussion  of  Results 

Results  obtained  in  the  four  problems  attempted  in 
this  chapter  are  analyzed  and  a  comparison  is  given  of 
the  two  basically  different  control  systems,  i.e.,  1)  an 
open  loop  control  system  and  2)  closed  loop  control  or 
specific  optimal  control  system. 

The  first  two  problems  analysed  represent  open  loop 
control.   Problem  1  is  simulated  in  Section  B.l  and  the 
operating  time  is  set  equal  to  one  second.   The  perform- 
ance criterion  used  here  could  be  interpreted  as  an  opti- 
mization of  control  and  movement,  limiting  erratic  flux 
changes  (10) . 

Another  useful  performance  index  is  the  subject  of 
problem  2  wherein  the  integral  of  the  error  squared  of 
the  reactor  power  and  its  desired  value  is  minimized  along 
with  the  integral  of  the  squared  of  the  reactivity  function. 
In  this  problem  the  control  function  required  to  achieve 
the  desired  state  for  example  n(t£)  =  an^  is  dependent 
on  the  value  of  a.   The  higher  the  value  of  a,  the  faster 
the  control  action  and  the  faster  the  power  rise.   This 
variation  is  exhibited  in  Figure  10.   Also  it  is  observed 
that  the  value  of  the  cost  function  increases  with  in- 
creasing value  of  a  (Figure  12).   This  is  to  be  expected 
since  the  higher  value  of  a  would  imply  desirability  of 


68 


driving  the  error  to  zero  at  a  faster  rate.   Naturally  more 
control  is  needed  to  accomplish  this. 

Problems  3  and  4  represent  closed  loop  or  specific 
optimal  control.   The  specific  forms  chosen  for  the 
controller  are  : 


p^  =  a  +  A(n=nj) 
and   p^   =  A(n-nd)  +  B   j  (n-nj)dt 


Th^  performance  index  used  is  the  same  as  in  Problem  2. 

The  results  of  problem  3,  point  out  a  strong 
dependence  of  the  peak  power  attained  to  the  desired 
value  of  the  power  level.   This  result  led  to  the  invest- 
igation of  Problem  4.   In  this  problem,  one  observes  that 
the  cost  function  increases  with  increase  in  the  desired 
peak  power  (Figure  17),  a  result  which  is  intuitively 
obvious. 

In  problem  4,  the  control  laws  are  determined  to 
achieve  four  different  power  levels  and  it  is  found  in 
each  case  that  the  closed  loop  or  specific  optimal  control 
suffers  in  performance  in  the  sense  that  the  cost  function 
is  higher  than  in  the  open  loop  case  (Figure  21).   If  the 
values  of  the  control  parameters  A  and  B  are  fixed  as  those 
obtained  for  the  case  N^  =  lOnQ  and  the  reactor  is  brought 
to  various  power  levels,  then  the  dotted  curve  in  Figure  22 
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Figure  21.  Cost  function  versus  peak  power  attained 
for  open  and  closed  loop  optimal  control 
(Nuclear  Reactor  System; 
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Figure  22.   Total  cost  function  versus  peak  power 

attained  for  open  loop  control  and  cost 
variation  for  fixed  and  optimally  varying 
values  of  control  paran eters  for  closed 
loop  control  (Nuclear  Reactor  System) 
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is  obtained  for  the  total  cost  function.   Note  here  that  th< 
cost  function  is  higher  in  every  case.   These  results  com- 
pare favorably  with  those  obtained  by  Eisenberg  and  Sage 
(21)  power  level  trajectories  and  the  corresponding  time 
variation  of  the  control  functions  are  shown  in  Figures  23 
and  24.   Figure  25  depicts  the  variation  of  the  cost  func- 
tion versus  the  peak  power  desired  for  varying  values  of 
Figure  26  shows  the  variation  of  the  cost  function  with 
variation  in  system  parameters  (neutron  lifetime).   It  is 
noted  here  that  the  closed^loop  control  is  better  than  the 
open  loop  control  for  large  variation  in  lifetime  about  its 
normal  value. 

It  may  be  pointed  out  at  this  stage  that  the  coeffi- 
cients in  the  control  law  are  highly  dependent  upon  the 
final  time  t^,  the  system  parameters,  and  the  initial  state 
of  the  system.   The  performance  of  the  system,  thus  is  a 
function  of  the  form  chosen  for  the  control  law.   The  form 
of  the  control  law  is  very  instrumental  in  determining  how 
well  the  performance  of  the  s.o.c.  compares  to  that  of  the 
true  optimal. 

Figure  18  shows  the  variation  of  the  parameters  A 
and  B  in  the  control  law,  as  a  function  of  the  peak 
power.   Since  this  variation  is  considerable,  adaptive 
control  is  visualized  as  a  means  of  further  improving  the 
performance.   If  adaptive  control  is  used,  the  example  at 
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Figure  23,   Open  loop  and  closed  loop  cptirnal 
trajectories  -  N(t)  versus  time 
(Nuclear  Reactor  System) 
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Figure  24.  Open  loop  and  closed  loop  optimal  control 
function  to  obtain  ten  tirres  the  initial 
power  level  in  a  Nuclear  Reactor  System 
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Figure  25.   Cost  function  versus  peak  power  desired 
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hand  could  be  that  of  the  so  called  performance  criterion 
sensing  or  extremum  adaptive  system.   An  overall  performance 
index  is  measured  and  an  attempt  is  made  to  find  a  control 
law  (by  trial  and  error)  which  optimizes  this  performance 
index.   The  learning  states  here  are  the  parameters  des- 
cribing the  control  law,  and  the  learning  process  consists 
of  the  trial  and  error  adjustment  of  these  parameters.   If 
the  learning  states  can  be  changed  rapidly  and  accurately 
the  adaptive  controller  may  be  able  to  handle  nonlinear 
control  systems  -  the  learning  states  will  represent  the 
best  instantaneous  linear  approximation  to  the  non-linear 
system.   This  approach  is  currently  under  study  (22) . 

The  gradient  technique  employed  yields  good  re- 
sults.  The  accuracy  of  the  results  can  be  considerably 
improved  at  the  expense  of  more  analog  and  logic  elements 
and  adjusting  values  of  K"s.   The  convergence  of  the  intera- 
tions  is  very  sensitive  to  the  values  of  the  K's.   Larger 
values  of  K's  result  in  excessive  overshoot  about  the  de- 
sired trajectory.   Smaller  values  result  in  an  increased 
number  of  cycles  to  converge. 

The  Quasilinearization  technique  is  next  employed 
to  compute  the  results  already  obtained  in  Problem  1. 
Figures  6  and  7  indicate  close  agreement  between  the 
results  obtained  by  the  two  techniques. 


CHAPTER  III 

OPTIMAL  CONTROL  OF  A  NUCLEAR  ROCKET  ENGINE 

Introduction  and  Model  Description 

The  nuclear  rocket  engine  to  be  considered  is 
similar  to  the  solid-core,  25,000  megawatt,  million 
pound  thrust  model  developed  by  Smith  and  Stenning  (23, 
24,25,26).    The  system  of  equations  describing  such  a 
model  are  based  on  a  lumped  parameter  analysis.   Such  an 
analysis  does  not  describe  the  system  corrpletely  and 
furthermore  the  analysis  is  affected  by  a  lack  of  agree= 
ment  among  investigators  concerning  the  temperature  re- 
activity effects.   Stenning  and  Smith  (27  )  contend  that 
temperature  reactivity  is  directly  proportional  to  the 
square  root  of  the  core  exit  stagnation  temperature,  while 
Mohler  and  Perry  (  2  )  hold  that  temperature  reactivity 
is  directly  proportional  to  this  temperature.   Weaver  (  12) 
observes,  that  both  contentions  could  give  fairly  accurate 
results  simply  by  choosing  appropriate  reactivity  coeffi^ 
cients.   Recent  rocket  engine  tests  (28  )  however  indicate 
that  temperature  reactivity  may  be  proportional  to  square 
of  this  temperature.   All  this  adds  up  to  the  fact  that  the 
nuclear  rocket  engines  possess  characteristics  which  need 
exhaustive  analysis.   Thus,  since  the  system  is  not  well 


77 


78 


understood,  a  closed=loop  control  is  desired.   The  maximum 
principle  however,  leads  to  an  open«=>loop  control.   The 
closed=loop  control  is  a  complex  function  of  the  state 
variables  and  due  to  the  analog  computer  limitations,  the 
treatment  is  restricted  to  open  loop  control  only. 

Basic  components  of  the  engine  are  depicted  in 
Figure  27  .   Liquid  hydrogen,  stored  in  a  low  pressure 
tank,  is  pum.ped  to  high  pressure,  and  is  used  to  regenera  = 
tively  cool  the  nozzle  and  reflectoro   The  propellant  gas 
is  further  heated  as  it  passes  through  the  reactor  core 
and  is  exhusted  through  the  nozzle„   A  portion  of  the  hot 
gas  is  bled  off  from  the  main  stream,  and  is  used  to  power 
the  turbo  pum.p ,  a  feature  giving  rise  to  the  name,  "Bleed 
turbine  system."   The  bleed  flow,  and  consequently  the  sys= 
tern  pressure  are  controlled  by  the  valve  in  front  of  the 
turbine. 

The  following  equations  describe  the  Smith-Stenning 
nuclear  rocket  engine  in  terms  of  the  time  dependent  nor= 
malized  quantities.   An  account  of  the  normalization  P^o  ^ 
cedure  is  given  in  references  (23)  and  (29)  and  therefore 
will  not  be  described  here. 


■i  M5 


B 


N'  =  P'N»  -  N«  +  C»  (3a) 


1  r 


A 


C»  =  N«  -  C«  (3„2) 
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N»  -  PHT)^^'^  (3.3) 


p«  =  P'V  (T«  ) 


1/2  .lElli 


P    '     .    ■     ''     '  (T  )l/2  (3^4) 

>=  Pc  ■"  a;  (T«)^^2  ■'a'  ~°  (3^^) 


where   N»  =  77-   .  C«  =  ~ 
Nd         Cd 


^   --fe 


and 


Equations  (3,1)  and  (3,2)  are  the  reactor  kinetics 
equations.   Equation  (3..3)  gives  the  time  rate  of  change 
of  maximum,  core  surface  temperature,  and  equation  (3,4) 
gives  the  time  rate  of  change  of  the  core  inlet  stagna- 
tion pressure.   Reactivity  is  given  by  equation  (3^5) 

Table  6  contains  important  constants  for  the  nuclear 
rocket  engine.   Note  that  the  delayed  neutron  fraction  is 
smaller  than  the  value  typically  associated  with  a  uranium- 
235  reactor.   This  reduced  value  was  used  by  Smith  and 
Stenning  to  account  for  diffusion  of  some  precursors  into 
coolant  channels  and  their  subsequent  expulsion  from  the 
system  before  emission  of  delayed  neutrons. 
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TABLE    6 

rcx:ket  engine  constants 


Design  reactor  power  25,000  MW 

Design  core  surface  temperature  5,000°R 

Design  core  inlet  stagnation  pressure  1,500  psi 

Thermal  time  constant  0.68  sec 

Pressure  time  constant  2.50  sec 

Control  poison  time  constant  2.00  sec 

^^alve  time  constant  0.20  sec 

Decay  constant  (one  group)  0.10  sec"^ 

Temperature  reactivity  coefficient  -10~""  0R~^ 

Hydrogen  reactivity  coefficient  0,3  °R/psi 

Neutron  lifetime  10~  sec 

Delayed  neutron  fraction  0.005 

Thrust,  F  1,500.000  lb 

Specific  impusle  Isp  640  sec 
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The  specifications  in  Table  6  significantly  exceed 
those  of  nuclear  engines  presently  under  development  and 
will  probably  be  first  achieved  through  clustering  of  small 
engines.   However,  since  the  above  system  has  received  con- 
siderable attention  in  the  open  literature,  it  makes  it  an 
attractive  system,  for  an  optimal  control  study^ 

It  may  be  pointed  out  here  that  the  reactivity 
equation  (3.5)  is  modified,  particularly  with  respect  to 
the  temperature  coefficient  of  reactivity.   The  (T")^^^ 
term  is  replaced  by  T\  as  in  the  Mohler=Perry  (  2  ) 
formulationo   In  view  of  this  modification,  equation  (3.5) 
is  re=written  as  : 

The  above  system  of  equations  with  the  constants 
of  Table  6  reduce  to: 

N'  -  50  (  P«N«-N«  +  C»)  (3.7) 

C»  =  0.1  (N«  -  CO  (3.8) 

f'  =  1.471  (N»  -  P«(T«)^'^^)    ^  (3.9) 

(T»)V2^  (3^10) 


2 

P'  =:  0.4  (P«  V»  (T«)^/2  -  IP;)^  ,-) 


p,  =p,^  -  T«  +  18  ^]  (3.11) 

The  system  is  in  an  idle  condition  for  T^O,  the 
idle  or  the  initial  condition  being  -.T^  =  1000°R  and  P^  = 
275  psi.   The  value  setting  is  fixed  at  0.915.   The  full 
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power  condition  is  T  =  5000°R  and  P  =  1500  psi.   It  is 
desired  to  optimize  system  performance  such  as  to  mini= 
mize  the  following  integral. 


tr 

J  =  1/2  i    p2(t)dt 
o 


with  T(tf)  =  T^  and    ip,  <  p^^^ 

This  performance  index  can  be  interpreted  as  tending  to 
limit  the  maximum  reactor  period  as  is  easily  seen  by 
considering  equation  (3.7). 

Since  the  effect  of  delayed  neutrons  is  small  in  the 
interval  of  operation,  N«/N«  =  50  P « .   This  indeed  is  the 
inverse  reactor  period 

By  using  the  statement  of  the  maximum  principle 

as  in  Chapter  II,  the  following  two  point  boundary  value 

problem  is  obtained. 

N»  =  50(  P  «M«  -  N»  +  C") 

C»  =  0.1(N»  -  CM 

T«  =  1.47  (N»  -  P»(T«)^'^^) 

(P^^2 


P«  =  0.4(P«V«(T")^/2  -  T^^flT^) 

P»  =   p»^  -  T«  +  18  !^ 

Pi  =-50  P'P^  +  50  P^  -  Oa  P2  -  1.46  P3 

P2  =-50  P^  +  0.1  P2 

P3  =  50  N'Pi  +  900  ^'^^'2'  ^   ^•'^^^  PMT")"  -   aP4 
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P=  -  SON'P, 
with  boundary  conditions: 
N«  (o)  z  0.0816 
C«(o)  =  0.0816 
T« (o)  =  0.2 
P« (o)  =  0.1S3 
Pl(tf )  =  0 
P2(tf )  =  0 
P^itf)   =  0 
T(tf)  =  aT«(o) 
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Solution  Procedure:   Since  the  useful  life  of  a  nuclear 

rocket  per  shot  is  approximately  2  to  15  minutes,  the 

nuclear  second  stage  must  begin  to  supply  full  thrust 

almost  immediately,  following  a  short  start-up  period. 

Thus  the  rocket  reactor  must  be  brought  up  to  full  power, 

of  hundreds  of  megawatts,  in  a  very  short  time.   In 

light  of  this  fact  the  optirrization  problem  was  solved 

for  various  ranges  of  desired  final  time,  4  $  t^  ^  11  sec, 

desired  final  state  0.4  ^T. ^1.0  and  variations  in  the 

a 

system  parameters,  2  t- ^  <  t  ^10  t    and  5  T^p  ^  T^f;   10  t^^ 
^^ariable  Start-up  Tjme:  4  $  tf  $:  11  sees. 

The  final  time  t^,  in  which  the  reactor  J.s  brought 
to  full  power  is  varied  between  4  sec  and  11  sec.   Three 
different  trajectories,  as  shown  in  Figure  28  are  used 
in  the  optimization  procedure.   The  corresponding  control 
poison  reactivity  and  the  pressure  responses  are  shown 
in  Figures  28  and  29.   These  are  obtained  by  forcing  the 
system  to  follow  the  desired  trajectory  in  an  optimum 
manner.   It  is  noted  from  the  figure,  that  shorter  the 
temperature  response  time,  the  faster  is  the  control  rod 
withdrawal  and  the  faster  is  the  pressure  rise  in  the 
system.   The  rod  has  to  be  re-inserted  to  maintain  the 
final  temperature  at  the  desired  level.   This  re-insertion 
is  very  gradual  for  longer  temperature  response  times. 
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''ariable  desired  temperature:  0.4  ^  j.   .;  1.0 

The  final  time  tf,  is  arbitrarily  fixed  to  10  sec 
and  the  system  is  driven  to  three  different  desired 
average  core  temperatures  from  the  idle  condition  of  0.2 
(1000°P).   In  each  case  it  is  desired  to  drive  the  system 
to  some  final  state  in  10  seconds  such  that 


1/2  ;' 


2 

P     dt     is  minimum 


o 

The  optimal  temperature  response  and  the  corresponding 
control  poison  reactivity  response  for  T^  =  0.4  (2000°F) 
is  given  in  Figure  30.   Figure  31  shows  the  optimum  pres- 
sure response  for  this  case.   Figures  32  and  33  show  the 
optimum  temperature,  control  poison  reactivity,  neutron 
density,  and  pressure  responses  for  T^  =   0,7  (3500OF)  and 
Td  =  1.0  (SOOOOF). 

With  the  optimal  trajectories  determined  and  the 
cost  function  J  computed  for  these  trajectories,  it  is  now 
possible  to  make  a  plot  of  the  cost  function  J  versus 
desired  average  core  temperature  (Figure  34).   The  cost 
to  drive  the  system  to  desired  core  temperatures  inr 
creases  slowly  at  first  and  rapidly  after  Td  =  0.4. 
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Figure  28.   Control  poison  reactivity  and  tetrperature 

time  plots,  tf  =  4  sees, 6  sees,  and  11  sees 

(Nuclear  Rocket  Engine  System) 
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Figure  29.   Pressure-time  plot,  tf  =  4  sees,  6  sees  and 
11  sees  (Nuclear  Rocket  Engine  System) 
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Figure  30.   Optimal  control  poison  reactivity  and  tempera- 
ture-time plots,  Tjj  =  0.4  (Nuclear  Rocket 
Engine  System) 
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Figure  32.   Optimal  control  poison  reactivity  and  tempera- 


ture-time plots 
Rocket  Engine  System) 


0.7  and  1.0  (Nuclear 
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Figure  33.   Cptirral  pressure  and  neutron  density  versus 

time,  Jd  =  0.7  and  1.0  (Nuclear  Rocket  Engine 
Systerr) 
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Figure  34.  Cost  function,  J  versus  desired  average 
core  terrperature  (Nuclear  Rocket  Engine 
System) 
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3 )   Influence  of  Parameter  v^ariation  on  System  Response 

The  following  four  cases  were  examined  with  respect 
to  the  normal  temperature  and  pressure  time  constants,  the 
normal  values  being,  y.   =  0.68  sec  and   T ^n  ~  ^'^  sec. 
The  reactivity  input  used  was  optimal  one  obtained  for 
the  case  of  7^   =  1.0. 

a)  Tp  =  2  Jpr^   ,  :^  =   t„ 

b)  yp  =  lO-Ypn   .Tt  =   tn 

c)  J  p  =  Tpn    .  Tt  =  5  Jtn 

d)  Tp  =  Tpn   .  U  =  10 


tn 


The  results  of  the  variation  of  the  pressure  time 
constant  from  the  normal  value  is  illustrated  in  Figure  35. 
As  the  pressure  time  constant  is  increased,  the  temperature 
response  becomes  quite  sluggish.   The  corresponding  pres- 
sure response  is  shown  in  Figure  36. 

If  the  pressure  time  constant  is  held  constant  and 
the  temperature  time  constant  is  increased,  the  system 
response  is  not  influenced  as  greatly  as  in  the  previous 
case  (Figure  37).   Humphries  (30)  has  reported  similar 
results  on  model  reference  adaptive  control  of  a  nuclear 
rocket  engine.   In  light  of  the  results  obtained  here  and 
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Humphries'  results,  it  appears  that  the  next  appropriate 
step  would  be  to  drive  the  system  with  a  control  whose 
specific  form  is  known.   The  controller  gains  will  have 
to  be  varied,  for  the  system  to  optimally  adapt  to  the 
desired  response.   If  this  is  not  satisfactory  an  optimal 
adaptive  controller  will  have  to  be  used.   Such  a  scheme 
is  currently  under  development  (  22) . 
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Figure  35, 


Control  poison  reactivity  versus  time  and  average 
core  terrperature  versus  time  for  three  cases  of 
parameter  variations  (Nuclear  Rocket  Engine  System) 
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Figure  36.   Pressure  versus  time  response  for  three  cases 
of  parameter  variations,  (Nuclear  Rocket 
Engine  Systerr) 
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Figure  37,   Pressure  versus  time  and  average  core  temper- 
ature versus  time  for  three  cases  of  para- 
meter variations  (Nuclear  Rocket  Engine  System) 


CHAPTER  IV 

SUMMARY  AND  CONCLUSIONS 

Summary  and  Results  :  Optimal  open-loop  control  and  the 
specific  optimal  control  of  nuclear  reactor  processes 
in  a  nuclear  reactor  system  and  open-loop  control  in  a 
nuclear  rocket  engine  system,  has  been  analysed  in  this 
study.   Optimal  control  theory  is  applied  to  the  problems 
and  the  analog  computer  is  used  for  simulation  and  in  some 
cases  optimization  by  the  gradient  technique  is  used.   In 
certain  examples,  quasilinearization  technique  is  used  for 
the  solution  of  the  TPBVP.   In  all  of  these,  the  converg- 
ence is  obtained  in  no  more  than  four  iterations. 

In  the  examples  considered  here,  for  nuclear  re- 
actor system,  the  closed-loop  control  suffers  in  perform- 
ance as  compared  to  the  open-loop  control  in  that  the  cost 
function  is  higher.   This  is  in  general  agreement  with  phy- 
sical insight  and  previous  results  (21),  The  coefficients  A 
and  B  in  the  control  law  are  highly  dependent  upon  the  final 
time  t£.  the  system  parameters,  and  the  initial  and 
desired  final  state.   Thus,  it  is  desirable  to  vary  the 
closed-loop  controller  parameters  as  a  function  of  neutron 
density  and  allowed  start-up  time,  as  well  as  system  para- 
meters, if  performance  near  optimum  is  to  be  achieved  for 
varying  start-up  time  and  final  neutron  density.   This 
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suggests  the  desirability  of  an  adaptive  type  control,  as 
a  means  of  further  improving  the  performance  for  such 
applications.   Such  a  scheme  is  currently  under  investi- 
gation by  Sage  and  Ellis  (22, 3i). 

The  open-loop  control  problem  for  the  nuclear 
rocket  engine  system  was  solved  in  three  different  parts: 
a)  variation  in  the  start-up  time,  b)  variation  in  the 
desired  average  core  temperature,  c)  influence  of  parameter 
variation  on  system  response. 

The  results  obtained  appear  quite  reasonable.   For 
example,  in  the  first  part  it  is  observed  that  shorter  the 
start-up  time,  the  faster  the  control  rod  withdrawal  and 
the  faster  the  pressure  rise.   Moreover,  in  order  to  main= 
tain  the  system  at  the  desired  final  state,  the  control 
rod  has  to  be  re-inserted  and  this  re-insertion  is  very 
gradual  for  longer  start-up  times  as  opposed  to  shorter 
start-up  times.   In  the  second  part  it  is  seen  that  the 
cost  function,  J,  increases,  if  it  is  desired  to  drive 
the  system  to  higher  and  higher  power  levels.   In  the 
third  part,  it  is  observed  that  if  the  pressure  and 
the  temperature  time  constants  were  increased  from  their 
normal  values ,  the  temperature  response  became  somewhat 
sluggish.   In  the  light  of  this  observation  and  Humphries 
(29)  results,  it  appears  that  the  next  appropriate  step 
would  be  to  drive  the  system  with  a  control  whose  specific 
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form  is  known.   The  controller  gains  will  have  to  be 
varied,  for  the  system  to  optimally  adapt  to  the  desired 
response.   Then  try  optimum  adaptive  or  adaptive  control. 

Since  it  is  possible  to  reduce  the  start-up  time 
considerably  using  optimal  control  theory  and  still  main- 
tain smooth  transition  to  full  power,  this  can -result  in 
large  amounts  of  propellant  savings  and  hence  reduction  in 
payload.   Does  computation  of  optimal  control  and  its  in- 
clusion in  driving  the  system  involve  additional  equipment 
to  render  it  impractical?  This  is  not  considered  to  be 
the  case  as  it  is  anticipated  that  computational  equipment 
will  be  available  in  the  total  payload.   Also  in  this  study, 
the  control  valve  opening  is  fixed  at  0.915  as  opposed  to 
the  maximum  of  1.33  to  3.0  considered  by  Smith  and  Stenning 
(23)  and  1.10  considered  by  Humphries  (29).   This  would  in 
turn  lower  the  turbine  size  because  of  lower  flows  and  thus 
allow  additional  reduction  in  weight. 

The  results  obtained  here  are  therefore  quite  en- 
couraging.  However,  this  control  concept  must  be  evaluated 
within  the  framework  of  the  total  body  of  knowledge  per- 
taining to  the  development  of  nuclear  rocketry.   How  these 
short  start-up  times  affect  the  attitude  and  guidance  con- 
trol problems  and  various  other  questions  pertaining  to  the 
equipment  and  their  costs  cannot  be  answered  here.   The 
results  are  sufficiently  encouraging  to  prompt  further 
investigation. 
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Suggestions  for  Further  Study 

A  logical  sequel  to  this  work  would  involve  the 
study  of  optimal  adaptive  control  of  the  two  systems  in 
order  to  improve  the  system  performance.   This  is  currently 
undertaken  by  Sage  and  Ellis.   Moreover,  the  control 
emphasis  should  be  shifted  from  the  control  poison  to  the 
hydrogen  propellant. 

An  optimal  control  of  more  complicated  reactor 
nodels  and  other  constraints  should  be  studied,  since  the 
optimal  process  is  only  as  good  as  the  constraints  and  the 
model  used.   In  particular  the  nuclear  reactor  system  with 
spatj.al  dependence  should  be  given  consideration  and  the 
reactor  system  should  be  modified  to  include  effect  of 
reflector. 
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APPENDIX  A 
ANALOG  SIMULATION 
Circuit  Symbolism 

The  analog  computer  wiring  diagrams  for  the  simula- 
tion of  the  four  control  problems  is  shown  in  Figures  41 
and  42.  .   The  explanation  of  the  circuit  diagrams  is 
given  in  Figures  38  and  39.    The  special  circuits 
used  are  shown  in  Figure  40  .   The  values  of  the  potentio- 
meter settings  (K)  are  contained  in  tables  rather  than  in 
the  computer  diagram.   Only  an  identification  number  is 
used  in  the  computer  diagram.   Diodes  are  used  in  the 
square  root  circuits  and  in  limiting  circuits. 

Input  and  feedback  resistances  for  summing  amplifiers 
and  input  resistances  for  integrating  amplifiers  are  ex- 
pressed in  megohms.   Values  of  either  0.1  or  1.0  megohms 
are  used.   Feedback  capacitances  of  integrating  amplifiers 
are  expressed  in  microfarads,  with  values  of  0.1  land  1.0 
microfarads  being  used  here.   An  identification  number  is 
assigned  to  each  amplifier  in  the  computer  circuit  diagram. 
Special  circuits  for  square-root  and  division  are  formed  by 
suitable  combinations  of  quarter-square  multipliers  and 
summing  amplifiers  (Figure  40  ).   With  the  symbolism  thus 
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Figure  38.   Symbols  used  in  analog  computer  diagrams 
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Figure  39,   Logic  symbols  used  in  corrputer  diagrams 
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defined,  it  is  now  possible  to  describe  the  various  siniu= 
lations ,  performed  in  the  course  of  this  investigation. 

Problem  Scaling  (Nuclear  Reactor  Systefm) 

Amplitude  and/or  time  scaling  is  usually  a  necessary 
part  of  analog  computer  work  to  contain  amplitudes  of  all 
variables  to  within  +  100  volts  and  time  to  within  measur- 
able quantities  on  the  equipment  used.   The  scaling  used  for 
the  four  problems  is  given  in  Table   "^  .   The  bars  above  the 
letters  denote  scaled  quantities.   With  this  scaling  the 
TPBVP  equations  for  the  four  control  problems  become: 

Problem  1 

N  =-100  N^P^  -  6.4  N  +  IOC 
t  =   0,064N  -  O.IC 
Pi  =  IOC  NP^^  +  6,4P^-6.4P2 
P2  =  -0.1  P^  +  0.1  P2 
Pc  =  0.1  N  P^ 
with  boundary  conditions: 

N(o)  =  0.05      ;    C(o)  =  0.032 
Pj^(tf)  *  0  and  'P2(tf )  -  0 

Problem  2 

N  =  lOON^P^  -  6.4  N  +  10  C 
C  =  0.064N  -  0,1  C 


Ill 


P^  =      lOON  P^^   +    6,4    P^    -    6.4    P2   +    a    (N-Nj) 
P2  =  -0.1   P^   +   0.1    P2 
Pc  =-0.1    N  Pi 

J  =   1/2       J       p^2   dt   +    5x10°^    a    J       O.l(N-Nd)    dt 
o  o 

J  =  J^  +  5x10^"^  a      J2' 

with   boundary   conditions    same   as    in  problem   1. 

Problem   3: 

N=0.1    aN+   AN    (N=Nj)    -   0.64      N  +   C 
t  =  0.0064   N  -    0.01   C 

Pi  =  -0.01(   a   +   A^)(N-Nj)    -    0,1    a   P^  -  2A    (NP^) 
+   A  N^   Pi   +    0.64   P^    -   0.64   P2   -   0.001    a   A 
P2  =  -0.01   P^   +   0.01   P2 

P3  =  -0.01   A    (N-N^)^-   Wi    (N-Nj)    -   0.001    a    (N-Nj) 
P4  =  -10""^   a   -   0.1    N  P^ 
Pc  =  10"^   a   +   0.01   A    (N-Nj),    A  is    negative 

Problem  4; 

Fi  =  AN    (N-Fid)    +0.1    BxN  -   0.64   N  +   C 

C  =  0.0064   N  -    0.01   C 

X  =  0.1    (N-Nd) 

P^   =   -0.01    (A^+a  )(N-N^)    =    0.2   A    (lOMPj^)    +    ANdPi 

-   0.1    BxP^    +   0.64pj    -    0.64P2   -   O.lP^   -    10'^  ABx 
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?2   =  -0.01  P^  +  0.01  P2 

P3  =-0.0l(0.01B^)x  -  0,1  BNP-^  -  AB  lO'^  (N-N^) 
P4  =  -0.01  A(N-N^)2  -  NPx(N=N^)  -lO"^  B  x(N-Nd) 
Pp^  =  =0.01  A(N-Nd)^  -0.01  (O.lB)x^  -  xNPj^ 
P^  =  0.01  A(N-Nd)  +  0.01  (OaB)x 
A  and  B  are  both  negative 
The  boundary  conditions  on  the  state  variables  are: 
N(o)  =  0.05  ,  C(o)  =  0.032 
x(o)  =  0;   P3  (t^)  =  P4(tf)  =  P5(tf)  =  0 
P4(o)  =  P4(o) =  0 

r^  ^  -   f^f         2 

J  =  1/2   J    p^2  dt  +  5x10°^  O      J  0.1(N-N  )  dt 
o  o 

J  =  J,  +  5x10°^  a    Jo° 


h  =     1/2   1     pM 


where   ^^ 

'o 


and  J2'  =   J   0.1(N-N.)  dt 


With  the  system  equations  scaled,  and  the  initial 
conditions  determined,  it  is  now  possible  to  proceed  with 
the  analog  computer  simulation.   The  analog  computer  cir- 
cuits for  simulating  the  above  problems  are  represented 
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in  Figures  38  and  39.   Potentiometer  settings  are  listed  in 
Tables  8,  9,  10  and  11. 

It  should  be  mentioned  here  that  great  care  and  in- 
genuity is  needed  to  arrange  the  analog  circuit,  so  that 
the  noise  is  not  amplified  at  any  particular  stage.   This 
is  particularly  brought  about  by  the  fact  that  in  the 
square  root  circuit,  the  output  is  multiplied  by  a  factor 
of  ten  and  in  the  division  circuit  the  output  is  multi- 
plied by  100  . 

TABLE  7 
SCALE  FACTORS  FOR  THE  FOUR  COrJTRCL  PROBLEMS 


Problem 

Scaled 

Que 

intities 

Variables 

Problem  1 

Problem 

T^ 

Problem  3 

-pi^BTi^r^ 

n 

ION 

ION 

ION 

ION 

C 

10% 

lO^C 

lO^C 

10% 

Pi 

10-^p^ 

10-^Pl 

io-4p^ 

10-'^P2 

P2 

10"^ P2 

10"^P2 

10"'^P2 

10-^P2 

0 

10-65 

lo-^a 

lo'^a 

P3 

P3 

10-^P3 

P4 

h 

P4 

,A 

10-3  A 

10  \ 

a 

lC-3i 

t 

r 

T 

lO-ir 

10"^  r 

P5 

h 

B 

10-^B 

X 

lOX 
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TABLE  8 
POT  E  NTT  IOMETER  SETTINGS  FOR  CONTROL  PROBLEMS  1  AND  2 
(NUCLEAR  REACTOR  SYSTEM) 

Potentio-  Parameter  Description  +100v  Problem  1  Problem  2 
meter 

1 

2 

3  Initial  Condition  on  N  +100 

4  Initial  Condition  on  C  +100  0.032 
5 
6 
7 
8 


9       Weighting  factor  for 
square  of  the  neutron 
density  error  term  i.e, 


10 


0,64 

0.64 

1.0 

1.0 

0.05 

0.05 

0.032 

0.032 

0.10 

0.10 

0.064 

0.064 

1.0 

1.0 

0„64 

0.64 

variable 

1.0 

1.0 

11  Initial  Condition  on  P^^  +100  variable  variable 

12  Initial  Condition  on  P2  "'"lOO  variable  variable 

13  0.64  0.64 

14  0.10  OolO 

15  0.10  0.10 

16  -100  N^(var.)  N^(vai:) 

17  O.ia 

18  K,  , 


varies 

varies 

from 

from 

0.01  to 

0.01  to 

0.04 

0.04 
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TABLE  8  (CONTINUED) 


Potentio-  Parameter  Description  +100v  Problem  1  Problem  2 

meter  


19  Kj2  varies    varies 

from  0.01  from  0.01 
to  0.04   to  0.04 

20  Koi  varies    varies 

^■^  from  0.01  from  0.01 

to  0.04    to  Oo04 

varies    varies 
from  0.01  from  0.01 
to  0.04   to  0.04 


21      K22 
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TABLE  9 

POTENTIOMETER  SETTINGS  FOR  CONTROL  PROBLEM  2 
(NUCLEAR  REACTOR  SYSTEM) 


Pot 

#16, 

Nd  = 

=  0.5 

Pot  #9 
0.1 

Pot  Hn 

PlO 

Pot  #12 
P20 

Pot  #16 
Nd 

2.0 

0.97 

0.13 

0.50 

2.0 

0.69 

0.21 

0.30 

2.0 

0.60 

0.24 

0.25 

2.0 

0.60 

0.41 

0.10 
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Simulation  by  Gradient  Technique 

The  algorithm  used  in  this  technique  is  given  by 
equation  (2.18),   The  problem  is  described  on  page  12 
and  the  analog  circuit  is  shown  in  Figure  41.   it  must 
be  borne  in  mind  that  the  integrators  must  all  be  started 
together;  thus  the  operate  (op)  inputs  are  all  joined  on 
the  logic  patchboard.   Each  time  the  op  input  receives 
a  logic  one,  a  solution  n(t)  is  generated.   By  systematic- 
ally changing  P^,    a  search  for  P^^*,  which  produces  the 
desired  optimal  control  u*  (reactivity  in  this  case)  and 
trajectory  n*  can  be  made.   The  computer  control  dia- 
gram for  iterative  computation  is  shown  in  Figure  43 
whereas  the  updating  circuit  is  in  Figure  ^'^   .   Before 
the  start  of  the  computation  the  Rep-op  is  off  and  the 
integrators  in  the  analog  circuit  are  in  the  I.C.  or 
reset  mode.   The  initial  conditions  P^"'"  and  P2"'~  a^e 
now  set.   The  first  run  begins  when  Rep-op  pushbutton 
is  depressed.   At  t  =  t^ ,  the  comparator  output  triggers 
the  memory  units  and  transfers  the  first  boundary  point 
to  their  output  (i.e  to  the  output  of  the  second  ampli- 
fier in  the  track-and-hold  system  of  the  updating  circuit 
(Figure  -^4  ).  since  the  outputs  are  held  until  t  =  tf  (or 
the  interval  set  by  the  pulser),  there  is  ample  time  for 
read-out  (on  the  brush  recorder  or  oscilloscope),  before 
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the  next  tun  begins.   This  is  a  great  advantage  of  the 

analog=digital  computer  because  it  can  show  the  human 

operator  as  to  actually  what  is  happening  via  an  on-line 

picture  on  the  oscilloscope.   The  iterations  can  be 

watched  step  by  step  and  the  source  of  trouble  can  be 

predicted  with  ease. 

Suppose  for  example,  it  is  desired  to  satisfy  the 

following  end  conditions:   Mlt^)  =:  10  n^  and  N(tf )  =  0. 

Then  the  p^"*"^  for  which  the  above  conditions  are  satisfied 
joS 

for  various  iteration  step  sizes  K^^,    K12.  ^21,  and  K22  is 
given  in  "^able  12,   It  is  clear  from  the  table  that  large 
values  of  K  result  In  large  overshoots  in  the  desired 
trajectory„   Sm.aller  values  however  result  in  increased 
number  of  iterations  to  converge^ 
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TABLE  12 


P^°S  AS  A  FUNCTION  OF  ITERATION  STEP  FOR  PROBLEM  1 

JO 


(NUCLEAR  REACTOR  SYSTEM) 


Run 


Iteration  step,  K^^  =  K22  =  ^21  =  K22  =  0.01 


^10      ^20   Pfo      ^00     ^  °^     Comments 
I Iterations 


1  0.2169  0.1762  0.42  0.1762  19  No  overshoot 

2  0.19  0.1762  0.50  0.1762  37  No  overshoot 

3  0.25  0„1762  0.50  0.1762  35  No  overshoot 

4  0.30  0.1762  0.53  0,1762  48  No  overshoot 

5  0.35     0,1762,   -     0.1762   131      Problem  term- 

inated 

6  0.40     0.1762   -     0.1762  very  large  Problem  term= 

inated 

7  0.45     0.1762   -     0.1762  very  large  Problem  term- 
„  ^  ors  inated 

8.  0.52     0.1762   -     0.1762  very  large  Problem  term= 
^  ^  inated 

9.  0.55     0.1762   -     0.1762  very  large  Very  close  to 

desired  tra- 
jectory 


Iteration  step,  K^^  =  K^2  =  ^2^  =  K22  =  0.02 

11  0.19     0.1762  0.46  0.1762    17  No  overshoot 

12  0.21     0.1762   0,505  0.1762    15  No  overshoot 

13  0.25     0.1762  0.51  0.1762    13  No  overfehoot 

14  0.30     0.1762   0.54  0.1762    19  No  overshoot 

15  0.35     0.1762   0.48  0,1762   :  75  Problem  term- 

inated 
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TABLE  12  (ContinuecJ) 


Run   Iteration  step,  K^^j^  =  K22  =  K21  =  K22  =  0.02 

#   — 

P{n      ^on       P(o      Pqo     ^  °f     Comments 
t_ ^^    ^" ^"   Iterations 

16  C.40     C.1762  0.393  0.1762  very  large  Problem  term- 

inated 

17  0.45     0.1762   0.444  0.1762  very  large  Problem  termi- 

nateci 

18  0.50     0.1762  0.50   0.1762  very  large  Problem  termi- 

nateci 

19  0.55     0.1762  0.55   0.1762  very  large  Problem  termi= 

nated 


Iter 

ation  step. 

^11  = 

K12  =  K21 

=  '<22  = 

0.03 

21 

0.19 

0.1762 

0,48 

0.1762 

12 

No  overshoot 

22 

0.21 

0.1762 

0.50 

0.1762 

10 

No  overshoot 

23 

0.25 

0.1762 

0,505 

0.1762 

8 

No  overshoot 

24 

0.30 

0.1762 

0.50 

0.1762 

10 

No  overshoot 

25 

0,35 

0.1762 

0.48 

0.1762 

40 

4%   overshoot 

26 

0.40 

0.1762 

0.394 

0.1762  very  large 

- 

Iteration  step, 

Kll  = 

K12  =  K21 

=  •<22  = 

0.04 

31 

0.19 

0.1762 

0.48 

0.1762 

7 

No  overshoot 

32 

0.21 

0.1762 

0,52 

0.1762 

6 

2%   overshoot 

33 

0.25 

0.1762 

0.47 

0^1762 

5 

No  overshoot 

34 

0.30 

0.1762 

0.49 

0.1762 

4 

No  overshoot 

35 

0.35 

0.1762 

0.60 

0.1762 

24 

A%   overshoot 

36 

0.40 

0.1762 

0.386 

0.1762  very  large 

- 
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TABLE  12  (Continued) 


Run 

# 


Iteration  step,  K^^  =  Ki2  =  K21  =  K22  =  0.01 


^10     '^20     ^10     ^20     ^  °^     Comments 


Al 

0.46 

0.37 

0.46 

A2 

0.46 

0.35 

0.46 

A3 

0.46 

0.32 

0.46 

A4 

0.46 

0,30 

0.46 

A5 

0.46 

0.27 

0.46 

A6 

0.46 

0.25 

0.46 

A7 

0.46 

0.20 

0.46 

A8 

0.46 

0.176 

0,46 

Iterations 

-    very  large 

"        65     No  overshoot 
(not  plotted) 
50     No  overshoot 
_  (not  plotted) 

50     No  overshoot 
_  (not  plotted) 

50     No  overshoot 
(not  plotted) 
0.16      44     20^  overshoot 

~  -  excessive 
overshoot 

~  -  excessive 
overshoot 


Iteration  step,  K^^  =  Ki2  =  K21  =  K22  =  0.02 


A12 

0.46 

0.35 

0.46 

- 

37 

No  overshoot 

A13 

0.46 

0.32 

0.46 

- 

33 

(not  plotted) 
No  overshoot 

A14 

0.46 

0.30 

0.46 

0.18 

28 

(not  plotted) 
No  overshoot 

A15 

0.46 

0.27 

0.46 

0.175 

35 

No  overshoot 

A16 

0.46 

0.25 

0.46 

0.155 

19 

30%  overshoot 

A17 

0.46 

0.20 

0.46 

0.13 

8 

excessive 

A18 

0.46 

0.176 

0.46 

- 

2 

overshoot 
excessive 
overshoot 
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TABLE  12  (Continued) 


Run   Iteration  step,  K^   =  K2^2  -   ^^21  =  ^^22  -  ^.03 

#    

pf^      pin     p!^    p!_     #  of     Comments 


10       20     '^10    '^20 


Iterations 


A21 

0.46 

0,37 

0.46 

0.16 

A22 

0.46 

0.35 

0.46 

0.14 

A23 

0.46 

0.32 

0.46 

0.16 

A24 

0.46 

0.30 

0.46 

0.16 

A25 

0.46 

0.27 

0.46 

0.16 

A26 

0.46 

0.25 

0.46 

0.16 

A27 

0.46 

0.20 

0.46 

- 

37 

14%  overshoot 

22 

38%  overshoot 

21 

12%  overshoot 

17 

8%  overshoot 

18 

8%  overshoot 

11 

20%  overshoot 

_ 

excessive 

overshoot 

Iteration  step,  K,,  =  K12  =  ^2^  =  K22 


0.04 


A30 

0.46 

0.37 

0.46 

0.17 

A31 

0.46 

0.35 

0.46 

0.18 

A32 

0.46 

0.32 

0.46 

0.165 

A33 

0.46 

0.30 

0.46 

0.16 

A34 

0.46 

0.25 

0.46 

- 

A35 

0.46 

0.20 

0.46 

- 

22 

8%  overshoot 

15 

2%  overshoot 

14 

14%  overshoot 

12 

20%  overshoot 

- 

excessive 

overshoot 

— 

excessive 

overshoot 

Iteration  step,  K^^^  =  K12  =  K21  =  K22  =  0.01 


0.46 

0.18 

0.50 

0.15 

8 

34%  overshoot 

0.40 

0.18 

0.44 

0.16 

7 

28%  overshoot 

0.35 

0.18 

0.39 

0,145 

7 

24%  overshoot 
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TABLE   12   (Continued) 


Run 

Itera 

ition  step, 

K 

11  = 

Ki2  = 

K21  =  K22  = 

0,01 

^10 

^^0 

pf 

UO 

^lo 

#  of 
Iterations 

Comments 

B4 

0o30 

0.18 

0 

.35 

0.14 

7 

2251^  overshoot 

B5 

0.25 

0.18 

0 

.31 

0.12 

8 

28%  overshoot 

B6 

0.20 

0.18 

0 

.28 

0.11 

13 

20%  overshoot 

B7 

0.20 

0.25 

- 

- 

- 

Does  not 
iterate 

B8 

0.25 

0.25 

~ 

■~ 

~ 

Does  not 
iterate 

B9 

0.30 

0.25 

0 

.41 

- 

19 

14%  overshoot 

BIO 

0.35 

0.25 

0 

.43 

0.155 

18 

12%  overshoot 

Bll 

0.40 

0.25 

0 

.495 

0,175 

30 

No  overshoot 

B12 

0.46 

0.25 

0 

.50 

0.18 

79 

8%  overshoot 

B13 

0.46 

0.30 

- 

- 

63 

No  overshoot 
(not  plotted) 

B14 

0.40 

0.30 

"~ 

43 

No  overshoot 
(not  plotted) 

B15 

0.35 

0.30 

"" 

^ 

35 

No  overshoot 
(not  plotted) 

B16 

0.30 

0.30 

" 

" 

~" 

Does  not 
iterate 

B17 

0.25 

0.30 

"" 

"~ 

~" 

Does  not 
iterate 

B18 

0.20 

0.30 

~" 

"~ 

~ 

Does  not 
iterate 

B19 

0.20 

0.35 

~ 

~ 

~" 

Does  not 
iterate 

B20 

0.25 

0.35 

^ 

"■ 

"" 

Does  not 
iterate 

B21 

0.30 

0.35 

~ 

-~ 

— 

Does  not 
iterate 

322 

0.35 

0.35 

"~ 

~ 

-~ 

Does  not 
iterate 

B23 

0.40 

0.35 

"" 

~~ 

— 

Does  not 
iterate 

B24 

0.46 

0.35 

" 

"" 

very  large 

No  overshoot 
(not  plotted) 
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TABLE  12  (Continued) 


Run 

# 


Iteration  step,  K^-^   =  Kx2  =  K21  =  K22  =  0-02 

^lo  "20        ^^0       4o     iJ.°^,o„3  """"-^"^^ 

Does  not 

iterate 

Does  not 

iterate 

12%  overshoot 

No  overshoot 

excessive 
overshoot 
no  overshoot 

A%   overshoot 

No  overshoot 

No  overshoot 

Does  not 

iterate 

Does  not 

iterate 

Does  not 

iterate 

Does  not 

iterate 

Does  not 

iterate 

Does  not 

iterate 

Does  not 

iterate 

Does  not 

iterate 

No  overshoot 

(not  plotted) 

Compare  with 

B24. 


C7 

0.20 

0.25 

- 

- 

— 

C8 

0.25 

0.25 

- 

- 

- 

C9 

0.30 

0.25 

0.39 

0.16 

9 

CIO 

0.35 

0,25 

- 

- 

8 

Cll 

0.40 

0.25 

- 

- 

12 

C12 

0.46 

0.25 

0.55 

0.18 

29 

CI  3 

0.46 

0.30 

0.59 

0.1625 

35 

C14 

0.40 

0.30 

0.53 

0.1775 

18 

C15 

0.35 

0.30 

0.4825 

0.16 

23 

C16 

0.30 

0.30 

- 

- 

- 

C17 

0.25 

0.30 

- 

- 

- 

C18 

0.20 

0.30 

- 

- 

- 

C19 

0.20 

0,35 

- 

- 

- 

C20 

0.25 

0.35 

- 

- 

- 

C21 

0.30 

0.35 

- 

- 

- 

C22 

0.35 

0.35 

- 

- 

- 

C23 

0.40 

0,35 

- 

- 

- 

C24 

0.46 

0.35 

- 

- 

52 
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Figure  43.   Computer  control  diagram 
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Figure  44.   Mpdating  circuit 
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APPENDIX  B 
COMPUTER  CODES 
1.   Solution  of  TPBVP  via  Quasilinearization 

The  FORTRAN  program  and  the  various  subroutines 
which  produced  some  of  the  results  in  Section  B.l  are 
given  in  the  pages  that  follow.   A  brief  summary  of  the 
code  is  given  below. 

After  the  necessary  data  have  been  read  in, the 
initial  approximation  is  generated  by  integrating  the 
non-linear  differential  equations.   The   (n+l)st  approxi- 
mation is  obtained  by  integrating  the  particular  and  homo- 
geneous equations,  determining  the  unknown  constants, 
and  forming  the  linear  combination  of  the  particular  and 
homogeneous  solutions.   This  (n+l)st  approximation  is 
printed  and  stored  as  the  nth  approximation  in  preparation 
for  the  next  iteration.   If  the  initial  values   of  the 
(n+l)st  approximation  agree  with  those  of  the  nth  approxi- 
mation to  within  five  or  more  decimal  places,  the  problem 
may  be  considered  solved. 

It  has  been  mentioned  before  that  the  convergence  is 
quadratic  ahd  hence  the  number  of  correct  digits  is  roughly 
doubled  with  each ,  iteration.   One  fact  that  stands  out 
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about  this  method  is  that  if  the  initial  approximation 
is  too  far  removed  from  the  solution,  the  iterations  may 
not  converge.   The  initial  approximation  used  in  this 
work  is  determined  from  the  analog  computer  results, 
■^his  is  found  to  be  very  reasonable  and  involves  large 
savings  in  computing  costs. 
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2.   Mimic  Code 

It  is  desirable  in  any  analog  computer  simulation  to 
have  a  means  for  checking  results.   A  digital  computer 
program  called  Mimic  (  32   )  is  used  for  this  purpose. 
Mimic,  a  digital  simulator  program   was  developed  for 
use  on  an  IBM  7090/7094  computer  with  a  FORTRAN  IV  IBSYS 
monitor.   The  University  of  Florida  Computing  Center 
personnel  modified  the  program  for  use  with  the  IBM  709 
computer. 

Mimic's  input  language  endows  the  digital  computer, 
with  the  parallel  nature  of  the  analog  computer  and  at  the 
same  time  eliminates  amplitude  and  tim.e  scaling.   If  the 
various  analog  operations  are  represented  by  black  boxes, 
then  the  equivalent  box-oriented  circuit  diagram  may  be 
used  for  writing  the  mimic  program.   The  following  example 
will  illustrate  this  fact.   Suppose  that  the  following 
equation  is  to  be  solved: 

X  +  X  +  X  =  0 
where  x(o)  =  2  and  x(o)  =  0 

To  solve  the  equation  on  the  analog  computer,  it  is  re- 
written in  the  following  form 

X  =  -X  -  X 
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Figure  45,     Box-oriented  mimic  diagram 
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The  box-oriented  circuit  diagram  is  given  in  Figure  45 
with  the  initial  conditions  specified  on  the  top  of  the 
boxes . 

Figure  46  illustrates  the  format  in  which  Mimic 
instructions  are  punched  on  an  IBM  card.   "Results"  are 
punched  beginning  in  column  10,  with  "expressions" 
beginning  in  column  19.   The  variable  names  2DX  and  IDX 
represent  the  second  and  first  derivatives  respectively, 
while  X  represents  the  desired  solution.   The  IMT  function 
integrates  the  first  argument  with  respect  to  the  inde- 
pendent variable,  and  the  initial  condition  is  the  second 
argument.   Mimic  employs  a  high  order  Runge-Kutta  sub- 
routine for  integration.   Note  that  the  Mimic  instructions 
proceed  very  much  like  the  corresponding  analog  computer 
program.   With  appropriate  control  instructions,  the  solu- 
tion X  may  be  obtained  as  a  function  of  the  independent 
variable. 

Mimic  contains  a  large  number  of  function  subroutines, 
including  logical  functions.   The  ease  with  which  its  in- 
structions can  be  presented  and  its  great  similarity  to 
analog  programming  make  it  an  ideal  instrument  for  checking 
analog  computer  scaling  and  setup.   Some  of  the  analog  com- 
puter results  checked  by  this  code  are  given  in  the  MIMIC- 
source  language  program  in  the  pages  that  fellow. 
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MIWIC   PROGRAM   NEI 


riRCN   FUTX   OPTIMISATION    (N  R   S) 
CON ( A, B. ALPHA) 
PAR(XND) 


SOLVE  EQHATIONS 

IDXN 

1000. *RHOC*XN-6.4*XN+0. 1*C 

XN 

INT (IDXN,. 5) 

IDC 

6.4*XN-  .1*C 

C 

INT(1DC,32.0) 

IDX 

XN=XND 

X 

INT (IDX,. 0) 

RHOC 

A*(XN=XND)+B*X 

HPOL 

RHOC*RH(X 

SPOL 

XN-XND 

\rpOL 

SPOL*SPOL 

BPCL 

. 5*HP0L+ . 5*ALPHA*VP0L 

XJ 

INT(BPOL,.0) 

FINISH  STATEMENT 

FIN(T,1.0) 

HEADER  STATEME\n- 

HDR(TIME) 

HDR 

HDR(XN,C.XJ) 

HDR 

HDR 

READOUT 

OUT(T) 

OUT 

OUT(XN,C,XJ) 

an 

OUT 

END 
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UIUIC   PROGRAM   NEirPRON   FL'JX   CPTIMISATION    (N   R    E) 


CON(^^DT,DTMIN) 
CFN{21.) 

XF 

SOLVE   EQUATIONS 

TP 

FIIN(XF,T) 

B 

TP 

A 

SQR(B) 

IDP 

.4♦P■»^A♦^^=.4*P*P/A 

P 

INT (IDP,. 183) 

TDOT 

der(b,t,o.o) 

TDOO 

TDOT 

XN 

TD00/1.47+P*A 

IDC 

0.1*XN-0,1*C 

C 

INT(1DC,.0816) 

XNDT 

DER(XN,T,0.0) 

XNDO 

x^roT 

RHO 

XND)/(50.^XN)+1.0-C/XN 

RH(X 

RH0+B-18.*P/B 

HPOL 

.5*RH0*RH0 

COST 

INT (HPOL, 0.0) 

FINISH   STATEMENT 

FIN  T,10,0) 

HEADER   STATEMENT 

HDR   TIME) 

HDR 

HDR  XN,TP,P,C) 

HDR   RHO, RHOC. COST) 

HDR 

HDR 

READOirr 

our(T) 

01  IT 

oirr(XN,Tp,p,c) 

01  JT(  RHO,  RHOC,  COST) 

OUT 

OUT 

END 
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